      SUBROUTINE WEIPLT(X,N,GAMMA)
C
C     PURPOSE--THIS SUBROUTINE GENERATES A WEIBULL
C              PROBABILITY PLOT
C              (WITH TAIL LENGTH PARAMETER VALUE = GAMMA).
C              THE PROTOTYPE WEIBULL DISTRIBUTION USED
C              HEREIN IS DEFINED FOR ALL POSITIVE X,
C              AND HAS THE PROBABILITY DENSITY FUNCTION
C              F(X) = GAMMA * (X**(GAMMA-1)) * EXP(-(X**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 WEIBULL 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  WEIBULL 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 WEIBULL 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, ALOG.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     REFERENCES--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 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--DECEMBER  1972. 
C     UPDATED         --MARCH     1975. 
C     UPDATED         --SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --FEBRUARY  1976. 
C
C---------------------------------------------------------------------
C
      DIMENSION X(1)
      DIMENSION Y(7500),W(7500)
      COMMON /BLOCK2/ WS(15000)
      EQUIVALENCE (Y(1),WS(1)),(W(1),WS(7501))
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 WEIPLT SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   17 FORMAT(1H , 98H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 WEIPLT SUBROUTINE IS OUTSIDE THE ALLOWABLE (1,,I6,16H) INTERVAL *
     1****)
   18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME
     1NT TO THE WEIPLT SUBROUTINE HAS THE VALUE 1 *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 WEIPLT 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
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     COMPUTE WEIBULL DISTRIBUTION ORDER STATISTIC MEDIANS
C
      DO100I=1,N
      W(I)=(-ALOG(1.0-W(I)))**(1.0/GAMMA)
  100 CONTINUE
C
C     PLOT THE ORDERED OBSERVATIONS VERSUS ORDER STATISTICS MEDIANS.
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
      CALL PLOT(Y,W,N)
      Q=.9975
      PP9975=(-ALOG(1.0-Q))**(1.0/GAMMA)
      Q=.0025
      PP0025=(-ALOG(1.0-Q))**(1.0/GAMMA)
      Q=.975
      PP975 =(-ALOG(1.0-Q))**(1.0/GAMMA)
      Q=.025
      PP025 =(-ALOG(1.0-Q))**(1.0/GAMMA)
      TAU=(PP9975-PP0025)/(PP975-PP025) 
      WRITE(IPR,105)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
      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
      SUM1=SUM1+(Y(I)-YBAR)*(Y(I)-YBAR) 
      SUM2=SUM2+(Y(I)-YBAR)*(W(I)-WBAR) 
      SUM3=SUM3+(W(I)-WBAR)*(W(I)-WBAR) 
  300 CONTINUE
      CC=SUM2/SQRT(SUM3*SUM1) 
      YSLOPE=SUM2/SUM3
      YINT=YBAR-YSLOPE*WBAR
      WRITE(IPR,305)CC,YINT,YSLOPE
C
  105 FORMAT(1H ,51HWEIBULL PROBABILITY PLOT WITH EXPONENT PARAMETER = ,
     1E17.10,1X,7H(TAU = ,E15.8,1H),11X,20HTHE SAMPLE SIZE N = ,I7)
  305 FORMAT(1H ,43HPROBABILITY PLOT CORRELATION COEFFICIENT = ,F8.5,5X,
     122HESTIMATED INTERCEPT = ,E15.8,3X,18HESTIMATED SLOPE = ,E15.8) 
C
      RETURN
      END 
      SUBROUTINE SORT(X,N,Y)
C
C     PURPOSE--THIS SUBROUTINE SORTS (IN ASCENDING ORDER)
C              THE N ELEMENTS OF THE SINGLE PRECISION VECTOR X
C              AND PUTS THE RESULTING N SORTED VALUES INTO THE
C              SINGLE PRECISION VECTOR Y.
C     INPUT  ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                                OBSERVATIONS TO BE SORTED. 
C                     --N      = THE INTEGER NUMBER OF OBSERVATIONS
C                                IN THE VECTOR X. 
C     OUTPUT ARGUMENTS--Y      = THE SINGLE PRECISION VECTOR
C                                INTO WHICH THE SORTED DATA VALUES
C                                FROM X WILL BE PLACED.
C     OUTPUT--THE SINGLE PRECISION VECTOR Y
C             CONTAINING THE SORTED
C             (IN ASCENDING ORDER) VALUES
C             OF THE SINGLE PRECISION VECTOR X.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--THE DIMENSIONS OF THE VECTORS IL AND IU 
C                   (DEFINED AND USED INTERNALLY WITHIN
C                   THIS SUBROUTINE) DICTATE THE MAXIMUM
C                   ALLOWABLE VALUE OF N FOR THIS SUBROUTINE.
C                   IF IL AND IU EACH HAVE DIMENSION K,
C                   THEN N MAY NOT EXCEED 2**(K+1) - 1.
C                   FOR THIS SUBROUTINE AS WRITTEN, THE DIMENSIONS
C                   OF IL AND IU HAVE BEEN SET TO 36,
C                   THUS THE MAXIMUM ALLOWABLE VALUE OF N IS
C                   APPROXIMATELY 137 BILLION.
C                   SINCE THIS EXCEEDS THE MAXIMUM ALLOWABLE
C                   VALUE FOR AN INTEGER VARIABLE IN MANY COMPUTERS,
C                   AND SINCE A SORT OF 137 BILLION ELEMENTS
C                   IS PRESENTLY IMPRACTICAL AND UNLIKELY,
C                   THEN THERE IS NO PRACTICAL RESTRICTION
C                   ON THE MAXIMUM VALUE OF N FOR THIS SUBROUTINE.
C                   (IN LIGHT OF THE ABOVE, NO CHECK OF THE 
C                   UPPER LIMIT OF N HAS BEEN INCORPORATED
C                   INTO THIS SUBROUTINE.)
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--THE SMALLEST ELEMENT OF THE VECTOR X
C              WILL BE PLACED IN THE FIRST POSITION
C              OF THE VECTOR Y,
C              THE SECOND SMALLEST ELEMENT IN THE VECTOR X
C              WILL BE PLACED IN THE SECOND POSITION
C              OF THE VECTOR Y, ETC.
C     COMMENT--THE INPUT VECTOR X REMAINS UNALTERED.
C     COMMENT--IF THE ANALYST DESIRES A SORT 'IN PLACE',
C              THIS IS DONE BY HAVING THE SAME
C              OUTPUT VECTOR AS INPUT VECTOR IN THE CALLING SEQUENCE. 
C              THUS, FOR EXAMPLE, THE CALLING SEQUENCE
C              CALL SORT(X,N,X)
C              IS ALLOWABLE AND WILL RESULT IN
C              THE DESIRED 'IN-PLACE' SORT.
C     COMMENT--THE SORTING ALGORTHM USED HEREIN
C              IS THE BINARY SORT.
C              THIS ALGORTHIM IS EXTREMELY FAST AS THE
C              FOLLOWING TIME TRIALS INDICATE.
C              THESE TIME TRIALS WERE CARRIED OUT ON THE
C              UNIVAC 1108 EXEC 8 SYSTEM AT NBS
C              IN AUGUST OF 1974.
C              BY WAY OF COMPARISON, THE TIME TRIAL VALUES
C              FOR THE EASY-TO-PROGRAM BUT EXTREMELY
C              INEFFICIENT BUBBLE SORT ALGORITHM HAVE
C              ALSO BEEN INCLUDED--
C              NUMBER OF RANDOM        BINARY SORT       BUBBLE SORT
C               NUMBERS SORTED
C                N = 10                 .002 SEC          .002 SEC
C                N = 100                .011 SEC          .045 SEC
C                N = 1000               .141 SEC         4.332 SEC
C                N = 3000               .476 SEC        37.683 SEC
C                N = 10000             1.887 SEC      NOT COMPUTED
C     REFERENCES--CACM MARCH 1969, PAGE 186 (BINARY SORT ALGORITHM
C                 BY RICHARD C. SINGLETON).
C               --CACM JANUARY 1970, PAGE 54.
C               --CACM OCTOBER 1970, PAGE 624.
C               --JACM JANUARY 1961, PAGE 41.
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         --NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DIMENSION X(1),Y(1)
      DIMENSION IU(36),IL(36) 
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)GOTO50
      IF(N.EQ.1)GOTO55
      HOLD=X(1)
      DO60I=2,N
      IF(X(I).NE.HOLD)GOTO90
   60 CONTINUE
      WRITE(IPR, 9)HOLD
      DO61I=1,N
      Y(I)=X(I)
   61 CONTINUE
      RETURN
   50 WRITE(IPR,15) 
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,18) 
      Y(1)=X(1)
      RETURN
   90 CONTINUE
    9 FORMAT(1H ,108H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT (A VECTOR) TO THE SORT   SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6
     1H *****)
   15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 SORT   SUBROUTINE IS NON-POSITIVE *****)
   18 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE SECOND INPUT ARGUME
     1NT TO THE SORT   SUBROUTINE HAS THE VALUE 1 *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
C     COPY THE VECTOR X INTO THE VECTOR Y
      DO100I=1,N
      Y(I)=X(I)
  100 CONTINUE
C
C     CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED
C
      NM1=N-1
      DO200I=1,NM1
      IP1=I+1
      IF(Y(I).LE.Y(IP1))GOTO200
      GOTO250
  200 CONTINUE
      RETURN
  250 M=1 
      I=1 
      J=N 
  305 IF(I.GE.J)GOTO370
  310 K=I 
      MID=(I+J)/2
      AMED=Y(MID)
      IF(Y(I).LE.AMED)GOTO320 
      Y(MID)=Y(I)
      Y(I)=AMED
      AMED=Y(MID)
  320 L=J 
      IF(Y(J).GE.AMED)GOTO340 
      Y(MID)=Y(J)
      Y(J)=AMED
      AMED=Y(MID)
      IF(Y(I).LE.AMED)GOTO340 
      Y(MID)=Y(I)
      Y(I)=AMED
      AMED=Y(MID)
      GOTO340
  330 Y(L)=Y(K)
      Y(K)=TT
  340 L=L-1
      IF(Y(L).GT.AMED)GOTO340 
      TT=Y(L)
  350 K=K+1
      IF(Y(K).LT.AMED)GOTO350 
      IF(K.LE.L)GOTO330
      LMI=L-I
      JMK=J-K
      IF(LMI.LE.JMK)GOTO360
      IL(M)=I
      IU(M)=L
      I=K 
      M=M+1
      GOTO380
  360 IL(M)=K
      IU(M)=J
      J=L 
      M=M+1
      GOTO380
  370 M=M-1
      IF(M.EQ.0)RETURN
      I=IL(M)
      J=IU(M)
  380 JMI=J-I
      IF(JMI.GE.11)GOTO310
      IF(I.EQ.1)GOTO305
      I=I-1
  390 I=I+1
      IF(I.EQ.J)GOTO370
      AMED=Y(I+1)
      IF(Y(I).LE.AMED)GOTO390 
      K=I 
  395 Y(K+1)=Y(K)
      K=K-1
      IF(AMED.LT.Y(K))GOTO395 
      Y(K+1)=AMED
      GOTO390
      END 
      SUBROUTINE UNIMED(N,X)
C
C     PURPOSE--THIS SUBROUTINE GENERATES THE N ORDER STATISTIC MEDIANS
C              FROM THE UNIFORM (RECTANGULAR)
C              DISTRIBUTION ON THE UNIT INTERVAL (0,1).
C              THIS DISTRIBUTION HAS MEAN = 0.5
C              AND STANDARD DEVIATION = SQRT(1/12) = 0.28867513.
C              THIS DISTRIBUTION HAS THE PROBABILITY
C              DENSITY FUNCTION F(X) = 1.
C              THIS SUBROUTINE IS A SUPPORT SUBROUTINE FOR
C              ALL OF THE PROBABILITY PLOT SUBROUTINES
C              IN DATAPAC; IT IS RARELY USED BY THE
C              DATA ANALYST DIRECTLY.
C              A PROBABILITY PLOT FOR A GENERAL DISTRIBUTION
C              IS A PLOT OF THE ORDERED OBSERVATIONS VERSUS 
C              THE ORDER STATISTIC MEDIANS FOR THAT DISTRIBUTION.
C              THE I-TH ORDER STATISTIC MEDIAN FOR A GENERAL
C              DISTRIBUTION IS OBTAINED BY TRANSFORMING
C              THE I-TH UNIFORM ORDER STATISTIC MEDIAN
C              BY THE PERCENT POINT FUNCTION OF THE DESIRED 
C              DISTRIBUTION--HENCE THE IMPORTANCE OF BEING ABLE TO
C              GENERATE UNIFORM ORDER STATISTIC MEDIANS.
C              IT IS OF THEROETICAL INTEREST TO NOTE THAT
C              THE I-TH UNIFORM ORDER STATISTIC MEDIAN
C              IN A SAMPLE OF SIZE N IS IDENTICALLY THE
C              MEDIAN OF THE BETA DISTRIBUTION
C              WITH PARAMETERS I AND N-I+1.
C     INPUT  ARGUMENTS--N      = THE DESIRED INTEGER NUMBER 
C                                OF UNIFORM ORDER STATISTIC MEDIANS
C                                TO BE GENERATED. 
C     OUTPUT ARGUMENTS--X      = A SINGLE PRECISION VECTOR
C                                (OF DIMENSION AT LEAST N)
C                                INTO WHICH THE GENERATED
C                                UNIFORM ORDER STATISTIC MEDIANS
C                                WILL BE PLACED.
C     OUTPUT--THE N ORDER STATISTIC MEDIANS
C             FROM THE RECTANGULAR DISTRIBUTION ON (0,1).
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     REFERENCES--FILLIBEN, 'THE PROBABILITY PLOT CORRELATION COEFFICIENT
C                 TEST FOR NORMALITY', TECHNOMETRICS, 1975, PAGES 111-117.
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         --SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DIMENSION X(1)
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)GOTO50
      IF(N.EQ.1)GOTO55
      GOTO90
   50 WRITE(IPR, 5) 
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR, 8) 
   90 CONTINUE
    5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 UNIMED SUBROUTINE IS NON-POSITIVE *****)
    8 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE FIRST  INPUT ARGUME
     1NT TO THE UNIMED 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 MEDIANS FOR THE FIRST AND LAST ORDER STATISTICS
C
      X(N)=0.5**(1.0/AN)
      X(1)=1.0-X(N) 
C
C     DETERMINE IF AN ODD OR EVEN SAMPLE SIZE
C
      NHALF=(N/2)+1 
      NEVODD=2*(N/2)
      IF(N.NE.NEVODD)X(NHALF)=0.5
      IF(N.LE.3)RETURN
C
C     COMPUTE THE MEDIANS FOR THE OTHER ORDER STATISTICS
C
      GAM=0.3175
      IMAX=N/2
      DO100I=2,IMAX 
      AI=I
      IREV=N-I+1
      X(I)=(AI-GAM)/(AN-2.0*GAM+1.0)
      X(IREV)=1.0-X(I)
  100 CONTINUE
C
      RETURN
      END 
      SUBROUTINE PLOT(Y,X,N)
C
C     PURPOSE--THIS SUBROUTINE YIELDS A ONE-PAGE PRINTER PLOT
C              OF Y(I) VERSUS X(I).
C     INPUT ARGUMENTS--Y      = THE SINGLE PRECISION VECTOR OF
C                               (UNSORTED OR SORTED) OBSERVATIONS
C                               TO BE PLOTTED VERTICALLY.
C                    --X      = THE SINGLE PRECISION VECTOR OF
C                               (UNSORTED OR SORTED) OBSERVATIONS
C                               TO BE PLOTTED HORIZONTALLY. 
C                    --N      = THE INTEGER NUMBER OF OBSERVATIONS
C                               IN THE VECTOR Y.
C     OUTPUT--A ONE-PAGE PRINTER PLOT OF Y(I) VERSUS X(I).
C     PRINTING--YES.
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--VALUES IN THE VERTICAL AXIS VECTOR (Y)
C              OR THE HORIZONTAL AXIS VECTOR (X) WHICH ARE
C              EQUAL TO OR IN EXCESS OF 10.0**10 WILL NOT BE
C              PLOTTED.
C              THIS CONVENTION GREATLY SIMPLIFIES THE PROBLEM
C              OF PLOTTING WHEN SOME ELEMENTS IN THE VECTOR Y
C              (OR X) ARE 'MISSING DATA', OR WHEN WE PURPOSELY
C              WANT TO IGNORE CERTAIN ELEMENTS IN THE VECTOR Y
C              (OR X) FOR PLOTTING PURPOSES (THAT IS, WE DO NOT
C              WANT CERTAIN ELEMENTS IN Y (OR X) TO BE PLOTTED).
C              TO CAUSE SPECIFIC ELEMENTS IN Y (OR X) TO BE 
C              IGNORED, WE REPLACE THE ELEMENTS BEFOREHAND
C              (BY, FOR EXAMPLE, USE OF THE   REPLAC   SUBROUTINE)
C              BY SOME LARGE VALUE (LIKE, SAY, 10.0**10) AND
C              THEY WILL SUBSEQUENTLY BE IGNORED IN THE PLOT
C              SUBROUTINE.
C     REFERENCES--NONE.
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         --OCTOBER   1974. 
C     UPDATED         --NOVEMBER  1974. 
C     UPDATED         --JANUARY   1975. 
C     UPDATED         --JULY      1975. 
C     UPDATED         --SEPTEMBER 1975. 
C     UPDATED         --OCTOBER   1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --FEBRUARY  1976. 
C     UPDATED         --FEBRUARY  1977. 
C
C---------------------------------------------------------------------
C
      CHARACTER*4 IGRAPH
      CHARACTER*4 SBNAM1,SBNAM2
      CHARACTER*4 ALPH11,ALPH12,ALPH21,ALPH22,ALPH31,ALPH32 
      CHARACTER*4 BLANK,HYPHEN,ALPHAI,ALPHAX
      CHARACTER*4 ALPHAM,ALPHAA,ALPHAD,ALPHAN,EQUAL
      DIMENSION Y(1)
      DIMENSION X(1)
      DIMENSION YLABLE(11)
      COMMON /BLOCK1/ IGRAPH(55,130)
C
      DATA SBNAM1,SBNAM2/'PLOT','    '/ 
      DATA ALPH11,ALPH12/'FIRS','T   '/ 
      DATA ALPH21,ALPH22/'SECO','ND  '/ 
      DATA ALPH31,ALPH32/'THIR','D   '/ 
      DATA BLANK,HYPHEN,ALPHAI,ALPHAX/' ','-','I','X'/
      DATA ALPHAM,ALPHAA,ALPHAD,ALPHAN,EQUAL/'M','A','D','N','='/
C
      IPR=6
      CUTOFF=(10.0**10)-1000.0
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      WRITE(IPR,998)
      IF(N.LT.1)GOTO52
      GOTO54
   52 WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,15)ALPH31,ALPH32,SBNAM1,SBNAM2
      WRITE(IPR,20)N
      WRITE(IPR,5)
      RETURN
   54 CONTINUE
      IF(N.EQ.1)GOTO56
      GOTO58
   56 WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,15)ALPH31,ALPH32,SBNAM1,SBNAM2
      WRITE(IPR,22)N
      WRITE(IPR,5)
      RETURN
   58 CONTINUE
C
      HOLD=Y(1)
      DO60I=2,N
      IF(Y(I).NE.HOLD)GOTO62
   60 CONTINUE
      WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,15)ALPH11,ALPH12,SBNAM1,SBNAM2
      WRITE(IPR,30)HOLD
      WRITE(IPR,5)
      RETURN
   62 CONTINUE
      HOLD=X(1)
      DO64I=2,N
      IF(X(I).NE.HOLD)GOTO66
   64 CONTINUE
      WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,15)ALPH21,ALPH22,SBNAM1,SBNAM2
      WRITE(IPR,30)HOLD
      WRITE(IPR,5)
      RETURN
   66 CONTINUE
C
      DO76I=1,N
      IF(Y(I).LT.CUTOFF)GOTO78
   76 CONTINUE
      WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,15)ALPH11,ALPH12,SBNAM1,SBNAM2
      WRITE(IPR,32) 
      WRITE(IPR,33)CUTOFF
      WRITE(IPR,5)
      RETURN
   78 CONTINUE
      DO80I=1,N
      IF(X(I).LT.CUTOFF)GOTO82
   80 CONTINUE
      WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,15)ALPH21,ALPH22,SBNAM1,SBNAM2
      WRITE(IPR,32) 
      WRITE(IPR,33)CUTOFF
      WRITE(IPR,5)
      RETURN
   82 CONTINUE
C
      N2=0
      DO96I=1,N
      IF(Y(I).LT.CUTOFF.AND.X(I).LT.CUTOFF)GOTO98 
      GOTO96
   98 N2=N2+1
      IF(N2.GE.2)GOTO99
   96 CONTINUE
      WRITE(IPR,5)
      WRITE(IPR,10) 
      WRITE(IPR,18)ALPH11,ALPH12,ALPH21,ALPH22
      WRITE(IPR,19)SBNAM1,SBNAM2
      WRITE(IPR,40) 
      WRITE(IPR,41)N2
      WRITE(IPR,5)
      RETURN
   99 CONTINUE
C
    5 FORMAT(1H ,'**************************************************',
     1'********************') 
   10 FORMAT(1H ,'                   FATAL ERROR                    ')
   15 FORMAT(1H ,'THE ',A4,A4,' INPUT ARGUMENT TO THE ',A4,A4,
     1' SUBROUTINE')
   18 FORMAT(1H ,'THE ',A4,A4,', AND ',A4,A4)
   19 FORMAT(1H ,'INPUT ARGUMENTS TO THE ',A4,A4,' SUBROUTINE')
   20 FORMAT(1H ,'IS NON-NEGATIVE (WITH VALUE = ',I8,1H))
   22 FORMAT(1H ,'HAS THE VALUE 1')
   30 FORMAT(1H ,'HAS ALL ELEMENTS = ',E15.8)
   32 FORMAT(1H ,'HAS ALL ELEMENTS IN EXCESS OF THE CUTOFF')
   33 FORMAT(1H ,'VALUE OF ',E15.8)
   40 FORMAT(1H ,'ARE SUCH THAT TOO MANY POINTS HAVE BEEN', 
     1' EXCLUDED FROM THE PLOT.')
   41 FORMAT(1H ,'ONLY ',I3,' POINTS ARE LEFT TO BE PLOTTED.')
C
C-----START POINT-----------------------------------------------------
C
C     DETERMINE THE VALUES TO BE LISTED ON THE LEFT VERTICAL AXIS
C
      DO200I=1,N
      IF(Y(I).GE.CUTOFF)GOTO200
      IF(X(I).GE.CUTOFF)GOTO200
      YMIN=Y(I)
      YMAX=Y(I)
      GOTO250
  200 CONTINUE
  250 DO300I=1,N
      IF(Y(I).GE.CUTOFF)GOTO300
      IF(X(I).GE.CUTOFF)GOTO300
      IF(Y(I).LT.YMIN)YMIN=Y(I)
      IF(Y(I).GT.YMAX)YMAX=Y(I)
  300 CONTINUE
      DO400I=1,9
      AIM1=I-1
      YLABLE(I)=YMAX-(AIM1/8.0)*(YMAX-YMIN)
  400 CONTINUE
C
C     DETERMINE THE VALUES TO BE LISTED ON THE BOTTOM HORIZONTAL AXIS 
C     DETERMINE XMIN, XMAX, XMID, X25 (=THE 25% POINT), AND 
C     X75 (=THE 75% POINT)
C
      DO600I=1,N
      IF(Y(I).GE.CUTOFF)GOTO600
      IF(X(I).GE.CUTOFF)GOTO600
      XMIN=X(I)
      XMAX=X(I)
      GOTO650
  600 CONTINUE
  650 DO700I=1,N
      IF(Y(I).GE.CUTOFF)GOTO700
      IF(X(I).GE.CUTOFF)GOTO700
      IF(X(I).LT.XMIN)XMIN=X(I)
      IF(X(I).GT.XMAX)XMAX=X(I)
  700 CONTINUE
      XMID=(XMIN+XMAX)/2.0
      X25=0.75*XMIN+0.25*XMAX 
      X75=0.25*XMIN+0.75*XMAX 
C
C     BLANK OUT THE GRAPH
C
      DO1100I=1,45
      DO1200J=1,109 
      IGRAPH(I,J)=BLANK
 1200 CONTINUE
 1100 CONTINUE
C
C     PRODUCE THE VERTICAL AXES
C
      DO1300I=3,43
      IGRAPH(I,5)=ALPHAI
      IGRAPH(I,109)=ALPHAI
 1300 CONTINUE
      DO1400I=3,43,5
      IGRAPH(I,5)=HYPHEN
      IGRAPH(I,109)=HYPHEN
 1400 CONTINUE
      IGRAPH(3,1)=EQUAL
      IGRAPH(3,2)=ALPHAM
      IGRAPH(3,3)=ALPHAA
      IGRAPH(3,4)=ALPHAX
      IGRAPH(23,1)=EQUAL
      IGRAPH(23,2)=ALPHAM
      IGRAPH(23,3)=ALPHAI
      IGRAPH(23,4)=ALPHAD
      IGRAPH(43,1)=EQUAL
      IGRAPH(43,2)=ALPHAM
      IGRAPH(43,3)=ALPHAI
      IGRAPH(43,4)=ALPHAN
C
C     PRODUCE THE HORIZONTAL AXES
C
      DO1500J=7,107 
      IGRAPH(1,J)=HYPHEN
      IGRAPH(45,J)=HYPHEN
 1500 CONTINUE
      DO1600J=7,107,25
      IGRAPH(1,J)=ALPHAI
      IGRAPH(45,J)=ALPHAI
 1600 CONTINUE
      DO1700J=20,107,25
      IGRAPH(1,J)=ALPHAI
      IGRAPH(45,J)=ALPHAI
 1700 CONTINUE
C
C     DETERMINE THE (X,Y) PLOT POSITIONS
C
      RATIOY=40.0/(YMAX-YMIN) 
      RATIOX=100.0/(XMAX-XMIN)
      DO1800I=1,N
      IF(Y(I).GE.CUTOFF)GOTO1800
      IF(X(I).GE.CUTOFF)GOTO1800
      MX=RATIOX*(X(I)-XMIN)+0.5
      MX=MX+7
      MY=RATIOY*(Y(I)-YMIN)+0.5
      MY=43-MY
      IGRAPH(MY,MX)=ALPHAX
 1800 CONTINUE
C
C     WRITE OUT THE GRAPH
C
      DO2100I=1,45
      IP2=I+2
      IFLAG=IP2-(IP2/5)*5
      K=IP2/5
      IF(IFLAG.NE.0)WRITE(IPR,2105)(IGRAPH(I,J),J=1,109)
      IF(IFLAG.EQ.0)WRITE(IPR,2106)YLABLE(K),(IGRAPH(I,J),J=1,109)
 2100 CONTINUE
      WRITE(IPR,2107)XMIN,X25,XMID,X75,XMAX
C
 2105 FORMAT(1H ,20X,109A1)
 2106 FORMAT(1H ,F20.7,109A1) 
 2107 FORMAT(1H ,14X,F20.7,5X,F20.7,5X,F20.7,5X,F20.7,1X,F20.7)
  998 FORMAT(1H1)
C
      RETURN
      END 
