SUBROUTINE WEIB(X,N)
C
C     PURPOSE--THIS SUBROUTINE PERFOMS A WEIBULL DISTRIBUTION ANALYSIS
C              ON THE DATA IN THE INPUT VECTOR X. 
C              THIS ANALYSIS CONSISTS OF DETERMINING THAT PARTICULAR
C              WEIBULL DISTRIBUTION
C              WHICH BEST FITS THE DATA SET.
C              THE GOODNESS OF FIT CRITERION IS THE MAXIMUM PROBABILITY
C              PLOT CORRELATION COEFFICIENT CRITERION.
C              AFTER THE BEST-FIT DISTRIBUTION IS DETERMINED,
C              ESTIMATES ARE COMPUTED AND PRINTED OUT FOR THE
C              LOCATION AND SCALE PARAMETERS.
C              TWO PROBABILITY PLOTS ARE ALSO PRINTED OUT-- 
C              THE BEST-FIT WEIBULL PROBABILITY PLOT
C              AND AN EXTREME VALUE TYPE 1 PROBABILITY PLOT 
C              (THIS IS DUE TO THE FACT THAT AS THE WEIBULL PARAMETER 
C              GAMMA APPROACHES INFINITY, THE WEIBULL DISTRIBUTION
C              APPROACHES THE EXTREME VALUE TYPE 1 DISTRIBUTION).
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     OUTPUT--4 PAGES OF AUTOMATIC PRINTOUT.
C     PRINTING--YES.
C     RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N
C                   FOR THIS SUBROUTINE IS 7500.
C     OTHER DATAPAC   SUBROUTINES NEEDED--SORT, UNIMED, WEIPLT,
C                                         EV1PLT, PLOT.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT AND ALOG.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     REFERENCE--FILLIBEN (1972), 'TECHNIQUES FOR TAIL LENGTH
C                ANALYSIS', PROCEEDINGS OF THE EIGHTEENTH
C                CONFERENCE ON THE DESIGN OF EXPERIMENTS IN 
C                ARMY RESEARCH AND TESTING, PAGES 425-450.
C              --FILLIBEN, 'THE PERCENT POINT FUNCTION',
C                UNPUBLISHED MANUSCRIPT.
C              --JOHNSON AND KOTZ (1970), CONTINUOUS UNIVARIATE
C                DISTRIBUTIONS-1, PAGES 250-271.
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--JUNE      1972. 
C     UPDATED         --AUGUST    1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --MAY       1976. 
C
C---------------------------------------------------------------------
C
      CHARACTER*4 IFLAG1
      CHARACTER*4 IFLAG2
      CHARACTER*4 IFLAG3
C
      CHARACTER*4 BLANK
      CHARACTER*4 ALPHAM
      CHARACTER*4 ALPHAA
      CHARACTER*4 ALPHAX
      CHARACTER*4 ALPHAI
      CHARACTER*4 ALPHAN
      CHARACTER*4 ALPHAF
      CHARACTER*4 ALPHAT
      CHARACTER*4 ALPHAY
      CHARACTER*4 ALPHAG
      CHARACTER*4 EQUAL
C
      DIMENSION W(3000)
      DIMENSION X(1)
      DIMENSION Y(7500),Z(7500)
      DIMENSION GAMTAB(50),CORR(50)
      DIMENSION YI(50),YS(50),T(50)
      DIMENSION IFLAG1(50),IFLAG2(50),IFLAG3(50)
      DIMENSION AINDEX(50)
      COMMON /BLOCK2/ WS(15000)
      EQUIVALENCE (Y(1),WS(1)),(Z(1),WS(7501))
      DATA BLANK,ALPHAM,ALPHAA,ALPHAX/' ','M','A','X'/
      DATA ALPHAI,ALPHAN,ALPHAF,ALPHAT,ALPHAY/'I','N','F','T','Y'/
      DATA ALPHAG,EQUAL/'G','='/
      DATA GAMTAB(1),GAMTAB(2),GAMTAB(3),GAMTAB(4),GAMTAB(5),
     1GAMTAB(6),GAMTAB(7),GAMTAB(8),GAMTAB(9),GAMTAB(10),
     1GAMTAB(11),GAMTAB(12),GAMTAB(13),GAMTAB(14),GAMTAB(15),
     1GAMTAB(16),GAMTAB(17),GAMTAB(18),GAMTAB(19),GAMTAB(20),
     1GAMTAB(21),GAMTAB(22),GAMTAB(23),GAMTAB(24),GAMTAB(25)
     1/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,
     113.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25./
      DATA GAMTAB(26),GAMTAB(27),GAMTAB(28),GAMTAB(29),GAMTAB(30),
     1GAMTAB(31),GAMTAB(32),GAMTAB(33),GAMTAB(34),GAMTAB(35),
     1GAMTAB(36),GAMTAB(37),GAMTAB(38),GAMTAB(39),GAMTAB(40),
     1GAMTAB(41),GAMTAB(42)
     1/30.,35.,40.,45.,50.,60.,70.,80.,90.,100.,150.,200.,250.,
     1350.,500.,750.,1000./
      DATA T(1),T(2),T(3),T(4),T(5),T(6),T(7),T(8),T(9),T(10),
     1T(11),T(12),T(13),T(14),T(15),T(16),T(17),T(18),T(19),T(20)
     1/1.63474,1.36116,1.34278,1.35854,1.37836,1.39657,1.41225,1.42557,
     1 1.43690,1.44660,
     1 1.45496,1.46223,1.46860,1.47422,1.47921,1.48368,1.48769,1.49132,
     1 1.49461,1.49761/
      DATA T(21),T(22),T(23),T(24),T(25),T(26),T(27),T(28),T(29),T(30),
     1T(31),T(32),T(33),T(34),T(35),T(36),T(37),T(38),T(39),T(40),
     1T(41),T(42),T(43)
     1/1.50036,1.50288,1.50521,1.50736,1.50935,1.51748,1.52344,1.52798,
     1 1.53157,1.53447,
     1 1.53888,1.54206,1.54447,1.54636,1.54788,1.55248,1.55480,1.55620,
     1 1.55781,1.55902,
     1 1.55997,1.56044,1.62391/
      DATA AINDEX(1),AINDEX(2),AINDEX(3),AINDEX(4),AINDEX(5),
     1AINDEX(6),AINDEX(7),AINDEX(8),AINDEX(9),AINDEX(10),
     1AINDEX(11),AINDEX(12),AINDEX(13),AINDEX(14),AINDEX(15),
     1AINDEX(16),AINDEX(17),AINDEX(18),AINDEX(19),AINDEX(20),
     1AINDEX(21),AINDEX(22),AINDEX(23),AINDEX(24),AINDEX(25)
     1/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,
     113.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25./
      DATA AINDEX(26),AINDEX(27),AINDEX(28),AINDEX(29),AINDEX(30),
     1AINDEX(31),AINDEX(32),AINDEX(33),AINDEX(34),AINDEX(35),
     1AINDEX(36),AINDEX(37),AINDEX(38),AINDEX(39),AINDEX(40),
     1AINDEX(41),AINDEX(42),AINDEX(43),AINDEX(44),AINDEX(45),
     1AINDEX(46),AINDEX(47),AINDEX(48),AINDEX(49),AINDEX(50)
     1/26.,27.,28.,29.,30.,31.,32.,33.,34.,35.,36.,37.,38., 
     139.,40.,41.,42.,43.,44.,45.,46.,47.,48.,49.,50./
C
      IPR=6
      IUPPER=7500
      NUMDIS=43
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1.OR.N.GT.IUPPER)GOTO50
      IF(N.EQ.1)GOTO55
      HOLD=X(1)
      DO60I=2,N
      IF(X(I).NE.HOLD)GOTO90
   60 CONTINUE
      WRITE(IPR, 9)HOLD
      RETURN
   50 WRITE(IPR,17)IUPPER
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,18) 
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE WEIB   SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 WEIB   SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL *
     1****)
   18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME
     1NT TO THE WEIB   SUBROUTINE HAS THE VALUE 1 *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      AN=N
C
C     COMPUTE THE SAMPLE MINIMUM AND SAMPLE MAXIMUM
C
      XMIN=X(1)
      XMAX=X(1)
      DO140I=2,N
      IF(X(I).LT.XMIN)XMIN=X(I)
      IF(X(I).GT.XMAX)XMAX=X(I)
  140 CONTINUE
C
C     COMPUTE THE PROB PLOT CORRELATION COEFFICIENTS FOR THE VARIOUS VALUES
C     OF GAMMA
C
      CALL SORT(X,N,Y)
      CALL UNIMED(N,Z)
C
      DO100IDIS=1,NUMDIS
      IF(IDIS.EQ.NUMDIS)GOTO150
      A=GAMTAB(IDIS)
      DO110I=1,N
      W(I)=(-ALOG(1.0-Z(I)))**(1.0/A)
  110 CONTINUE
      GOTO170
  150 DO160I=1,N
      W(I)=-ALOG(ALOG(1.0/Z(I)))
  160 CONTINUE
C
  170 SUM1=0.0
      SUM2=0.0
      DO200I=1,N
      SUM1=SUM1+Y(I)
      SUM2=SUM2+W(I)
  200 CONTINUE
      YBAR=SUM1/AN
      WBAR=SUM2/AN
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      DO300I=1,N
      SUM2=SUM2+(Y(I)-YBAR)*(W(I)-WBAR) 
      SUM1=SUM1+(Y(I)-YBAR)*(Y(I)-YBAR) 
      SUM3=SUM3+(W(I)-WBAR)*(W(I)-WBAR) 
  300 CONTINUE
      SY=SQRT(SUM1/(AN-1.0))
      CC=SUM2/SQRT(SUM3*SUM1) 
      YSLOPE=SUM2/SUM3
      YINT=YBAR-YSLOPE*WBAR
      CORR(IDIS)=CC 
      YI(IDIS)=YINT 
      YS(IDIS)=YSLOPE
  100 CONTINUE
C
C     DETERMINE THAT DISTRIBUTION WITH THE MAX PROB PLOT CORR COEFFICIENT
C
      IDISMX=1
      CORRMX=CORR(1)
      DO400IDIS=1,NUMDIS
      IF(CORR(IDIS).GT.CORRMX)IDISMX=IDIS
      IF(CORR(IDIS).GT.CORRMX)CORRMX=CORR(IDIS)
  400 CONTINUE
      DO500IDIS=1,NUMDIS
      IFLAG1(IDIS)=BLANK
      IFLAG2(IDIS)=BLANK
      IFLAG3(IDIS)=BLANK
      IF(IDIS.EQ.IDISMX)GOTO550
      GOTO500
  550 IFLAG1(IDIS)=ALPHAM
      IFLAG2(IDIS)=ALPHAA
      IFLAG3(IDIS)=ALPHAX
  500 CONTINUE
C
C     WRITE OUT THE TABLE OF PROB PLOT CORR COEFFICIENTS FOR VARIOUS GAMMA
C
      WRITE(IPR,998)
      WRITE(IPR,305)
      WRITE(IPR,999)
      WRITE(IPR,310)N
      WRITE(IPR,311)YBAR
      WRITE(IPR,312)SY
      WRITE(IPR,313)XMIN
      WRITE(IPR,314)XMAX
      WRITE(IPR,999)
      WRITE(IPR,323)
      WRITE(IPR,324)
      WRITE(IPR,325)
      WRITE(IPR,999)
C
      NUMDM1=NUMDIS-1
      IF(NUMDM1.LT.1)GOTO850
      DO800I=1,NUMDM1
      WRITE(IPR,805)GAMTAB(I),CORR(I),IFLAG1(I),IFLAG2(I),IFLAG3(I),
     1YI(I),YS(I),T(I)
  800 CONTINUE
  850 I=NUMDIS
      WRITE(IPR,806)ALPHAI,ALPHAN,ALPHAF,ALPHAI,ALPHAN,ALPHAI,
     1ALPHAT,ALPHAY,CORR(I),IFLAG1(I),IFLAG2(I),IFLAG3(I),
     1YI(I),YS(I),T(I)
C
C     PLOT THE PROB PLOT CORR COEFFICIENT VERSUS GAMMA VALUE INDEX
C
      CALL PLOT(CORR,AINDEX,NUMDIS)
      WRITE(IPR,810)ALPHAG,ALPHAA,ALPHAM,ALPHAM,ALPHAA,EQUAL,
     1GAMTAB(1),GAMTAB(12),GAMTAB(23),GAMTAB(34), 
     1ALPHAI,ALPHAN,ALPHAF,ALPHAI,ALPHAN,ALPHAI,ALPHAT,ALPHAY
      WRITE(IPR,999)
      WRITE(IPR,812)
      WRITE(IPR,813)
C
C     IF THE OPTIMAL GAMMA IS FINITE, PLOT OUT THE WEIBULL
C     PROBABILITY PLOT FOR THE OPTIMAL VALUE
C     OF GAMMA.
C
      IF(IDISMX.LT.NUMDIS)CALL WEIPLT(X,N,GAMTAB(IDISMX))
C
C     PLOT OUT AN EXTREM VALUE TYPE 1 PROBABILITY PLOT
C     (WHICH IS IDENTICALLY A WEIBULL PROBABILITY 
C     WITH GAMMA = INFINITY)
C
      CALL EV1PLT(X,N)
C
  998 FORMAT(1H1)
  999 FORMAT(1H )
  305 FORMAT(1H ,40X,16HWEIBULL ANALYSIS)
  310 FORMAT(1H ,37X,20HTHE SAMPLE SIZE N = ,I7)
  311 FORMAT(1H ,34X,18HTHE SAMPLE MEAN = ,F14.7) 
  312 FORMAT(1H ,28X,32HTHE SAMPLE STANDARD DEVIATION = ,F14.7)
  313 FORMAT(1H ,32X,21HTHE SAMPLE MINIMUM = ,F14.7)
  314 FORMAT(1H ,32X,21HTHE SAMPLE MAXIMUM = ,F14.7)
  323 FORMAT(1H ,85H       WEIBULL          PROBABILITY PLOT     LOCATIO
     1N         SCALE       TAIL LENGTH)
  324 FORMAT(1H ,83H     TAIL LENGTH          CORRELATION        ESTIMAT
     1E        ESTIMATE       MEASURE)
  325 FORMAT(1H ,37H   PARAMETER (GAMMA)      COEFFICIENT)
  805 FORMAT(1H ,3X,F10.2,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5) 
  806 FORMAT(1H ,5X,8A1,13X,F8.5,1X,3A1,2X,F14.7,2X,F14.7,3X,F10.5)
  810 FORMAT(1H ,12X,5A1,1X,A1,F14.7,11X,F14.7,11X,F14.7,11X,F14.7,
     115X,8A1)
  812 FORMAT(1H ,96HTHE ABOVE IS A PLOT OF THE 46 PROBABILITY PLOT CORRE
     1LATION COEFFICIENTS (FROM THE PREVIOUS PAGE))
  813 FORMAT(1H ,16X,35HVERSUS THE 46 WEIBULL DISTRIBUTIONS)
C
      RETURN
      END