C     DPLOD2.F  - STAND ALONE FORTRAN CODE FOR LIMITS OF DETECTION
C                 METHOD USED BY THE ASTM E2677 STANDARD.
C
C     THESE CODES WERE EXTRACTED FROM THE DATAPLOT SOFTWARE.  THEY HAVE
C     BEEN MODIFIED SLIGHTLY (PRIMARILY THE I/O) FOR THIS STAND ALONE
C     CODE.  THE SPECIAL FUNCTION ROUTINES ARE ORIGINALLY FROM THE
C     SLATEC LIBRARY (THE TPPF ROUTINE CALLS DBETAI WHICH BRINGS IN AN
C     ADDITIONAL SEVEN SPECIAL FUNCTION ROUTINES). 
C
C     NOTE THAT WITH gfortran UNDER Linux, YOU CAN COMPILE THIS
C     ROUTINE WITH:
C
C        gfortran -c -D LINUX  dplod2.F
C
C     THIS WILL CREATE THE OBJECT FILE dplod2.o WHICH CAN BE LINKED
C     WITH YOUR DRIVER ROUTINE (THIS SHOULD ALSO WORK FOR Mac OSX
C     AS LONG AS YOU HAVE INSTALLED THE gfortran COMPILER).
C
C     FOR WINDOWS PLATFORMS, THE SYNTAX WILL DEPEND ON THE SPECIFIC
C     COMPILER USED.  HOWEVER, YOU WILL NEED TO USE THE PRE-PROCESSOR
C     DIRECTIVE
C
C            -D WIN32
C
C     IF YOU DO NOT WANT TO USE THE PRE-PROCESSOR DIRECTIVE, YOU CAN
C     COMMENT/UNCOMMENT THE APPROPRIATE "LINUX" OR "WIN32" LINES IN
C     THE INITMC ROUTINE BELOW.  NO OTHER SOURCE CODE CHANGES SHOULD BE
C     NEEDED.
C
C     THIS FILE CONTAINS THE FOLLOWING ROUTINES:
C
C     DPLOD2   - THIS IS THE ROUTINE THAT YOU WILL CALL
C     DPLOD3   - THIS ROUTINE IS CALLED BY DPLOD2 WHEN THE
C                "CLEAN" CASE IS DETECTED.
C
C     THE FOLLOWING UTILITLY ROUTINES USED BY DPLOD2, DPLOD3 ARE
C     ALSO INCLUDED.
C
C         1. CHSPFF    - PERCENT POINT FUNCTION OF CHI-SQUARE DISTRIBUTION
C         2. D9LGMC    - CORRECTION FACTOR FOR LOG GAMMA (USED BY TPPF)
C         3. DBETAI    - INCOMPLETE BETA FUNCTION (USED BY TPPF)
C         4. DCSEVL    - EVALUATE CHEBYCHEV SERIES (USED BY TPPF)
C         5. DLBETA    - NATURAL LOGARITHM OF THE COMPLETE BETA FUNCTION
C                        (USED BY TPPF)
C         6. DLNGAM    - LOG GAMMA FUNCTION (USED BY TPPF)
C         7. DLNREL    - EVALUATE LN(1+X) (USED BY TPPF)
C         8. DGAMMA    - GAMMA FUNCTION (USED BY TPPF)
C         9. DGAMLM    - USED BY DGAMMA (USED BY TPPF)
C        10. DISTIN    - RETURN DISTINCT ELEMENTS OF A VARIABLE
C        11. INITDS    - INITIALIZE CHEBYCHEV SERIES (USED BY TPPF)
C        12. INITMC    - DEFINE MACHINE CONSTANTS (FOR WINDOWS AND LINUX)
C        13. MEAN      - COMPUTE THE MEAN
C        14. NORPPF    - PERCENT POINT FUNCTION OF NORMAL DISTRIBUTION
C        15. PERCEN    - COMPUTE A SAMPLE PERCENTILE
C        16. SD        - COMPUTE THE STANDARD DEVIATION
C        17. SORT      - QUICK SORT
C        18. SORTC     - QUICK SORT AND CARRY ANOTHER VARIABLE
C        19. TPPF      - PERCENT POINT FUNCTION OF T DISTRIBUTION
C        20. VAR       - COMPUTE THE VARIANCE
C
      SUBROUTINE DPLOD2(Y,X,N,MAXOBV,IPR,
     1                  XIDTEM,TEMP1,TEMP2,TEMP3,YT,YM,
     1                  D,BA,BL,SU,T,P,LO,LU,S,UN,TO,UA,
     1                  SI,B1,YTT,YVAR,
     1                  NS,MU,NN,
     1                  ALPHA,BETA,GAMMA,PA,THRPER,CVUSER,
     1                  ILODCV,ILODOF,ILODST,ILODTA,ILODPC,
     1                  ALOD,ALODSE,CV,ALODLU,
     1                  IDEBUG,IERROR)
C
C     PURPOSE--PERFORM A LIMITS OF DETECTION ANALYIS.  SPECIFICALLY,
C              WE ARE CHECKING FOR A LOWER LIMIT OF DETECTION.
C              IT IS ASSUMED THAT WE CANNOT MEASURE THE MASS BELOW
C              A CERTAIN DETECTION LIMIT AND WE NEED TO DETERMINE
C              WHAT THAT DETECTION LIMIT IS.  THERE ARE 2 CASES:
C
C                  1) FOR A ZERO MASS, THE MEASUREMENT WILL ALWAYS
C                     RETURN A ZERO VALUE.  THIS WILL BE REFERRED
C                     TO AS THE "CLEAN" CASE".
C
C                  2) FOR A ZERO MASS, THE MEASUREMENT MAY BE A
C                     NON-ZERO VALUE.  THIS WILL BE REFERRED TO AS
C                     THE DIRTY CASE.
C
C              THE DATA CONSISTS OF A MASS VALUE ("X") AND A
C              MEASURED VALUE ("Y").  THE COMMAND SYNTAX IS:
C
C                  CLEAN LIMIT OF DETECTION  Y   X
C                  DIRTY LIMIT OF DETECTION  Y   X
C                  LIMIT OF DETECTION  Y   X
C
C              IF THE CLEAN/DIRTY CASE IS NOT SPECIFIED, ASSUME
C              THE DIRTY CASE.
C
C     SUBROUTINES NEEDED--NORPPF, CHSPPF
C     WRITTEN BY--ALAN HECKERT
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2899
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--2018/08
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES--------------
C
      CHARACTER*4 ILODCV
      CHARACTER*4 ILODOF
      CHARACTER*4 ILODST
      CHARACTER*4 ILODTA
      CHARACTER*4 ILODPC
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
C----------------------------------------------------------------
C
      REAL Y(*)
      REAL X(*)
      REAL TEMP1(*)
      REAL TEMP2(*)
      REAL TEMP3(*)
      REAL YT(*)
      REAL YM(*)
      REAL XIDTEM(*)
C
      REAL D(*)
      REAL BA(*)
      REAL BL(*)
      REAL SU(*)
      REAL T(*)
      REAL P(*)
      REAL LO(*)
      REAL LU(*)
      REAL S(*)
      REAL UN(*)
      REAL TO(*)
      REAL UA(*)
      REAL SI(*)
      REAL B1(*)
      REAL YTT(*)
      REAL YVAR(*)
C
      INTEGER NS(*)
      INTEGER MU(*)
      INTEGER NN(*)
C
      REAL ZERO
      REAL EPS
C
      DOUBLE PRECISION DSUM
C
      INTEGER NCHAR(8)
C
      CHARACTER*4 IQUAME
      CHARACTER*4 ISTEPN
      CHARACTER*11 IVALC
      CHARACTER*60 IFORMT
      CHARACTER*132 IHEAD
C
      LOGICAL IFLAGY
      LOGICAL IFLAGZ
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA ZERO /0.0/
      DATA EPS  /0.1E-7/
C
C-----START POINT------------------------------------------------
C
C     SET OUTPUT UNIT (6 FOR SCREEN, IF NOT SET TO 6, THEN OPEN
C     THE OUTPUT FILE TO "LOD.OUT").
C
CCCCC IPR=6
      IF(IPR.LT.6 .OR. IPR.GE.100)IPR=6
      IF(IPR.NE.6)THEN
        OPEN(IPR,FILE='LOD.OUT')
      ENDIF
C
      IERROR='NO'
      NUMDIG=5
C
C     SPECIFY METHOD FOR COMPUTING PERCENTILES (SEE COMMENTS
C     IN THE "PERCEN" ROUTINE FOR THE DEFINITIONS).  FOR LOD
C     PURPOSES, SHOULD NOT REALLY MAKE MUCH DIFFERENCE WHICH
C     METHOD YOU PICK.
C
      IQUAME='R6'
CCCCC IQUAME='R7'
CCCCC IQUAME='R9'
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF DPLOD2--')
        WRITE(IPR,52)N,ALPHA,BETA,GAMMA,PA
   52   FORMAT('N,ALPHA,BETA,GAMMA,PA = ',I8,4G15.7)
        WRITE(IPR,53)THRPER,CVUSER
   53   FORMAT('THRPER,CVUSER = ',2G15.7)
        DO55I=1,N
          WRITE(IPR,56)I,Y(I),X(I)
   56     FORMAT('I,Y(I),X(I) = ',I8,2G15.7)
   55   CONTINUE
      ENDIF
C
      CALL INITMC(IDEBUG,IPR)
      CPUMIN=-R1MACH(2)
      CPUMAX=R1MACH(2)
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               **  1) NEGATIVE VALUES FOR THE RESPONSE   **
C               **     ARE CONSIDERED ERRORS.             **
C               **  2) IF THERE IS NO "ZERO" MASS VALUE,  **
C               **     CONSIDER THIS AN ERROR FOR NOW.    **
C               ********************************************
C
      ISTEPN='1.1'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
  991 FORMAT('AT STEP ',A4)
C
      IF(N.LT.2)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
  101   FORMAT('***** ERROR IN LIMITS OF DETECTION--')
        WRITE(IPR,102)
  102   FORMAT('      THE NUMBER OF OBSERVATIONS FOR THE LIMITS OF ',
     1         'DETECTION')
        WRITE(IPR,103)N
  103   FORMAT('      IS LESS THAN 2;  THE ENTERED NUMBER OF ',
     1         'OBSERVATIONS = ',I8)
        WRITE(IPR,999)
        IERROR='YES'
        GOTO9000
      ENDIF
C
      YMIN=CPUMAX
      DO200I=1,N
        IF(ABS(X(I)).LE.EPS)X(I)=ZERO
        YMIN=MIN(YMIN,Y(I))
        IF(X(I).LT.ZERO)THEN
          WRITE(IPR,999)
          WRITE(IPR,101)
          WRITE(IPR,212)
  212     FORMAT('      A NEGATIVE VALUE WAS ENCOUNTERED FOR ',
     1           'THE MASS VALUE')
          WRITE(IPR,204)I,X(I)
  204     FORMAT('      ROW ',I8,' HAS THE VALUE ',G15.7)
          IERROR='YES'
          GOTO9000
        ENDIF
  200 CONTINUE
C
      IF(YMIN.LT.0.0)THEN
        YCONST=ABS(YMIN) + 1.0
        DO208I=1,N
          Y(I)=Y(I) + YCONST
  208   CONTINUE
      ELSE
        YCONST=0.0
      ENDIF
C
C               ********************************************
C               **  STEP 1B-                              **
C               **  CHECK FOR:                            **
C               **  1) ARE THERE ANY MEASUREMENTS FOR     **
C               **     MASS = 0?  IF NOT, THIS IS AN      **
C               **     ERROR (AT LEAST FOR NOW).          **
C               **  2) ARE ALL MEASUREMENTS FOR MASS = 0  **
C               **     EQUAL TO ZERO?  IF SO, USE THE     **
C               **     "CLEAN" ALGORITHM.                 **
C               **     IF THERE ARE LESS THAN 3 NON-ZERO  **
C               **     MEASUREMENTS FOR MASS = 0, THEN    **
C               **     REMOVE THESE FROM THE DATA AND USE **
C               **     THE CLEAN ALGORITHM.               **
C               ********************************************
C
C     MAY 2009: MICHAEL VERKOUTEREN HAS SUGGESTED THE FOLLOWING
C               MINIMUM REQUIREMENTS:
C
C               1) 10 OR MORE BLANK REPITITIONS
C               2) AT LEAST 3 LEVELS
C               3) AT LEAST 10 REPITITIONS FOR EACH LEVEL
C
C     JUNE 2009: MICHAEL VERKOUTEREN HAS MADE THE FOLLOWING ADDITIONAL
C                SUGGESTIONS:
C
C                1) IF > 90% OF MEASUREMENTS AT HIGHTEST LEVEL
C                   ARE ZERO, THEN REPORT ERROR.
C                2) IF DATA IN ALL LEVELS IS EITHER ALL ZERO'S OR
C                   ALL NON-ZEROS, THEN RECOMMEND MEASUREMENTS
C                   AT MASS VALUE NEAR MID-POINT TWO SUCCESSIVE
C                   MASS LEVELS WHERE LOWER IS ALL ZERO MEASUREMENTS
C                   AND UPPER IS ALL NON-ZERO MEASUREMENTS.
C
C                   09/2010: FOR CASE 2), MICHAEL WANTS THE
C                            RECOMMENDATION PRINTED, BUT HE STILL
C                            WANTS THE RESULTS REPORTED.
C
C     AUGUST 2009: MICHAEL VERKOUTEREN ADDED ADDITIONAL CRITIERION.
C
C                1) IF > 10% OF MEASUREMENTS ARE LESS THAN THE
C                   CRITICAL VALUE, THEN WE NEED TO REPORT AN
C                   ERROR THAT DATA NEEDS TO BE TAKEN AT A HIGHER
C                   LEVEL.
C
      ISTEPN='1.2'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
C
 1000 CONTINUE
      CALL DISTIN(X,N,XIDTEM,NDIST,IPR,IDEBUG,IERROR)
      CALL SORT(XIDTEM,NDIST,XIDTEM,IPR)
C
      IF(XIDTEM(1).NE.ZERO)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
        WRITE(IPR,1002)
 1002   FORMAT('      THE LIMIT OF DETECTION COMMAND REQUIRES')
        WRITE(IPR,1003)
 1003   FORMAT('      SOME MEASUREMENTS FOR MASS = 0.')
        WRITE(IPR,1004)XIDTEM(1)
 1004   FORMAT('      THE MINIMUM MASS VALUE  = ',G15.7)
        IERROR='YES'
        GOTO9000
      ELSEIF(NDIST.LT.3)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
        WRITE(IPR,1012)
 1012   FORMAT('      AT LEAST THREE DISTINCT LEVELS REQUIRED.')
        WRITE(IPR,1013)NDIST
 1013   FORMAT('      THE NUMBER OF LEVELS DETECTED WAS ',I5)
        IERROR='YES'
        GOTO9000
      ELSEIF(NDIST.EQ.N)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
        WRITE(IPR,1014)
 1014   FORMAT('      THE NUMBER OF DISTINCT MASS VALUES IS EQUAL TO ',
     1         'THE SAMPLE SIZE')
        IERROR='YES'
        GOTO9000
      ENDIF
C
C     MAY 2009: MOVE CODE TO DETERMINE THE NUMBER OF REPITITIONS
C               FOR EACH LEVEL HERE.  WE NEED THIS TO DETERMINE
C               IF WE HAVE SUFFICIENT REPITITIONS FOR EACH LEVEL
C               AND TO DETERMINE WHETHER WE NEED THE CLEAN OR
C               DIRTY CASE.
C
      ISTEPN='1.3'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
C
      CALL SORTC(X,Y,N,X,TEMP1,IPR)
      DO1018I=1,N
        Y(I)=TEMP1(I)
 1018 CONTINUE
C
C     MAY 2009: CHECK FOR MINIMUM NON-ZERO VALUE. DETERMINE
C               IF THIS IS AT A ZERO MASS CASE OR NOT.  THIS
C               IS USED TO DECIDE HOW TO COMPUTE H IN CODE
C               BELOW.
C
      YMIN=CPUMAX
      YMINMS=CPUMAX
      DO1022I=1,N
        IF(Y(I).GT.0.0)THEN
          IF(Y(I).LT.YMIN)THEN
            YMIN=Y(I)
            YMINMS=X(I)
          ENDIF
        ENDIF
 1022 CONTINUE
C
C     JULY 2009: DETERMINE "SECOND SMALLEST VALUE" FROM ALL THE
C                DATA, NOT JUST THE MASS = 0 DATA.
C
      YMIN2=CPUMAX
      DO1024I=1,N
        IF(Y(I).GT.YMIN)THEN
          IF(Y(I).LT.YMIN2)THEN
            YMIN2=Y(I)
          ENDIF
        ENDIF
 1024 CONTINUE
C
      ISTEPN='1.4'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
C
      DO1030I=1,NDIST
C
        HOLD=XIDTEM(I)
        NREP=0
        NZERO=0
        NPLUS=0
C
        DO1040J=1,N
          IF(X(J).EQ.HOLD)THEN
            NREP=NREP+1
            IF(Y(J).LE.0.0)THEN
              NZERO=NZERO+1
            ELSE
              NPLUS=NPLUS+1
            ENDIF
          ENDIF
 1040   CONTINUE
C
        IF(IDEBUG.EQ.'ON')THEN
          WRITE(IPR,1041)I,XIDTEM(I),NREP,NZERO,NPLUS
 1041     FORMAT('I,XIDTEM(I),NREP,NZERO,NPLUS=',I5,G15.7,3I5)
        ENDIF
C
        IF(NREP.LT.10)THEN
          WRITE(IPR,999)
          WRITE(IPR,101)
          WRITE(IPR,1032)
 1032     FORMAT('      AT LEAST 10 REPITITIONS ARE REQUIRED FOR ',
     1           'EACH LEVEL.')
          WRITE(IPR,1033)I,HOLD,NREP
 1033     FORMAT('      LEVEL ',I5,' (= ',G15.7,') ONLY HAS ',I5,
     1           'REPITITIONS.')
          IERROR='YES'
          GOTO9000
        ELSEIF(I.EQ.1)THEN
          NUMZER=NZERO
          NUMPOS=NPLUS
        ENDIF
C
 1030 CONTINUE
C
C     THE THRESHOLD FOR USING THE "CLEAN" ALGORITHM IS THAT
C     >= 90% OF THE BLANK CASES MUST BE ZERO.
C
C     SET AN ADDITIONAL THRESHOLD OF >= 50% OF BLANKS ARE
C     NON-ZERO.
C
C     SET IFLAG = 0 FOR CLEAN CASE, IFLAG = 1 FOR DIRTY CASE,
C     AND IFLAG = 2 FOR CASE WHERE WE NEED TO CHECK FOR SUFFICIENT
C     NUMBER OF BLANKS (REQUIRE AT LEAST 20).
C
      ATEMP=REAL(NUMZER+NUMPOS)*0.9 + 0.000001
      MINZER=INT(ATEMP)
      ATEMP=REAL(NUMZER+NUMPOS)*0.5 + 0.000001
      MINPOS=INT(ATEMP)
      IF(NUMZER.GE.MINZER)THEN
        IFLAG=0
      ELSEIF(NUMPOS.GE.MINPOS)THEN
        IFLAG=1
      ELSE
        IFLAG=2
      ENDIF
C
C     CHECK FOR SUFFICIENT NUMBER OF BLANK CASES
C
      IF(IFLAG.EQ.2 .AND. NUMZER+NUMPOS.LT.20)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
        WRITE(IPR,1052)
 1052   FORMAT('      WHEN THE PERCENTAGE OF NON-ZERO MEASUREMENTS ',
     1           'FOR BLANKS')
        WRITE(IPR,1053)
 1053   FORMAT('      IS LESS THAN 50%, AT LEAST 20 BLANK ',
     1         'MEASUREMENTS ARE REQUIRED.')
        NADD=20 - (NUMZER+NUMPOS)
        WRITE(IPR,1054)NADD
 1054   FORMAT('      AT LEAST ',I5,' ADDITIONAL BLANK MEASUREMENTS ',
     1         'ARE REQUIRED.')
        IERROR='YES'
        GOTO9000
      ENDIF
C
C     FOR CLEAN CASE, REMOVE ANY POSITIVE MEASUREMENTS
C
      IF(IFLAG.EQ.0 .AND.NUMPOS.GT.0)THEN
        ICNT=0
        DO1065I=1,N
          IF(X(I).EQ.XIDTEM(1) .AND. Y(I).GT.ZERO)THEN
            GOTO1065
          ELSE
            ICNT=ICNT+1
            X(ICNT)=X(I)
            Y(ICNT)=Y(I)
          ENDIF
 1065   CONTINUE
        N=ICNT
      ENDIF
C
      IF(IFLAG.EQ.0)THEN
        CALL DPLOD3(Y,X,N,MAXOBV,IPR,
     1              XIDTEM,TEMP1,TEMP2,TEMP3,YT,YM,
     1              D,BA,BL,SU,T,P,LO,LU,S,UN,TO,UA,SI,B1,YTT,YVAR,
     1              NS,MU,NN,
     1              ALPHA,BETA,GAMMA,PA,THRPER,CVUSER,
     1              ALOD,ALODSE,CV,ALODLU,
     1              YMINMS,YMIN,YMIN2,
     1              ILODCV,ILODOF,ILODST,ILODTA,ILODPC,
     1              IDEBUG,IERROR)
        GOTO9000
      ENDIF
C
      IF(THRPER.EQ.CPUMIN)THRPER=0.9
C
C               ***********************************************
C               **  STEP 2.1--                               **
C               **  PERFORM THE BASIC CALCULATIONS.          **
C               ***********************************************
C
      ISTEPN='2.1'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
C
C
C     ALPHA   = PROBABILITY OF FALSE POSITIVE READINGS
C               (NO SIGNAL, BUT ALARM SOUNDS)
C     BETA    = PROBABILITY OF FALSE NEGATIVE READINGS
C               (SIGNAL, BUT NO ALARM)
C     GAMMA   = DEFINES CONFIDENCE LIMIT FOR LOD
C               (ALSO DEFINES THE CONFIDENCE FOR THE TOLERANCE LIMIT)
C     PA      = DEFINES THE COVERAGE FACTOR FOR THE TOLERANCE LIMIT
C
      CALL NORPPF(1.0-ALPHA,ZA,IPR)
      CALL NORPPF(1.0-BETA,ZB,IPR)
      CALL NORPPF(1.0-PA,ZP,IPR)
      GA=GAMMA
C
      IF(IDEBUG.EQ.'ON')THEN
        ISTEPN='2.2'
        WRITE(IPR,991)ISTEPN
        WRITE(IPR,999)
        WRITE(IPR,61)
   61   FORMAT('***** FROM DPLOD2--')
        WRITE(IPR,62)ZA,AZ,ZB,ZP
   62   FORMAT('ZA,AZ,ZB,ZP = ',4G15.7)
      ENDIF
C
C     COMPUTE:
C
C       1)  NUMBER OF BLANK (ZERO) MASS VALUES           (M)
C       2)  NUMBER OF DISTINCT MASS VALUES               (NDIST)
C       3)  SAMPLE SIZES FOR EACH GROUP                  (NS)
C       4)  NUMBER OF ZERO RESPONSES IN EACH GROUP       (MU)
C       5)  NUMBER OF NON-ZERO RESPONSES IN EACH GROUP   (NN)
C       6)  MEAN OF NON-ZERO RESPONSES IN EACH GROUP     (YM)
C
C     FOR THE MASS = 0 CASE, COMPUTE SOME STATISTICS
C
      NFINAL=0
      IFLAGY=.FALSE.
      IFLAGZ=.TRUE.
      IALLZ=-99
      IALLP=-99
C
      DO1100I=1,NDIST
C
        D(I)=0.0
        BL(I)=0.0
        SU(I)=0.0
        T(I)=0.0
        P(I)=0.0
        IF(I.GT.1)THEN
          LO(I)=0.0
          LU(I)=0.0
          BA(I)=0.0
          UN(I)=0.0
          S(I)=0.0
          TO(I)=0.0
        ELSE
          LO(I)=99.0
          LU(I)=99.0
          BA(I)=99.0
          UN(I)=99.0
          S(I)=99.0
          TO(I)=99.0
        ENDIF
C
        NZERO=0
        NPLUS=0
        NABOVE=0
        NBELOW=0
        HOLD=XIDTEM(I)
        DSUM=0.0D0
        DO1110J=1,N
          IF(X(J).EQ.HOLD)THEN
            IF(Y(J).LE.0.0)THEN
              NZERO=NZERO+1
            ELSE
              NPLUS=NPLUS+1
              TEMP1(NPLUS)=Y(J)
              DSUM=DSUM + DBLE(Y(J))
              IF(I.GT.1 .AND. TEMP1(NPLUS).GE.CV)NABOVE=NABOVE+1
            ENDIF
C
            IF(I.EQ.NDIST)THEN
              IF(Y(J).LT.CV)NBELOW=NBELOW+1
            ENDIF
C
          ENDIF
C
 1110   CONTINUE
        NTOTZ=NZERO+NPLUS
        MU(I)=NZERO
        NN(I)=NPLUS
        NS(I)=NTOTZ
C
        IF(NZERO.GE.1 .AND. NPLUS.GE.1)IFLAGZ=.FALSE.
        IF(IFLAGZ)THEN
          IF(NZERO.EQ.NTOTZ)IALLZ=I
          IF(NPLUS.EQ.NTOTZ .AND. IALLP.LT.0)IALLP=I
        ENDIF
        IF(I.EQ.NDIST)THEN
          NTEMPZ=INT(0.90*REAL(NTOTZ)+0.1)
          NTEMPB=INT(0.10*REAL(NTOTZ)+0.1)
          IF(NZERO.GT.NTEMPZ)THEN
            IFLAGY=.TRUE.
            WRITE(IPR,999)
            WRITE(IPR,101)
            WRITE(IPR,7071)
 7071       FORMAT('      More than 90% of the measurements at the')
            WRITE(IPR,7072)
 7072       FORMAT('      highest mass level are zero.')
            WRITE(IPR,7073)
 7073       FORMAT('      The LOD analysis will not be performed.')
            WRITE(IPR,999)
            IERROR='YES'
            GOTO9000
          ELSEIF(NBELOW.GT.NTEMPB)THEN
            IFLAGY=.TRUE.
            WRITE(IPR,999)
            WRITE(IPR,101)
            WRITE(IPR,7076)
 7076       FORMAT('      More than 10% of the measurements at the')
            WRITE(IPR,7077)
 7077       FORMAT('      highest mass level are less than the ',
     1             'critical value..')
            WRITE(IPR,7073)
            WRITE(IPR,7078)
 7078       FORMAT('      It is recommended that you take ',
     1             'measurements at a higher mass level.')
            WRITE(IPR,999)
            IERROR='YES'
            GOTO9000
          ENDIF
        ELSEIF(I.EQ.2)THEN
          AVAL=REAL(NABOVE)/REAL(NTOTZ)
          IF(AVAL.GT.0.90)THEN
            WRITE(IPR,999)
            WRITE(IPR,7180)
 7180       FORMAT('***** WARNING IN LIMITS OF DETECTION--')
            WRITE(IPR,7181)
 7181       FORMAT('      More than 90% of the measurements at the ',
     1             'lowest non-zero')
            WRITE(IPR,7183)
 7183       FORMAT('      mass level are greater than the critical ',
     1             'value.')
            WRITE(IPR,7185)
 7185       FORMAT('      Although the LOD analysis will still be ',
     1             'performed,')
            WRITE(IPR,7187)
 7187       FORMAT('      it is recommended that measurements be ',
     1             'taken')
            AJUNK=XIDTEM(2)/2.0
            WRITE(IPR,7189)AJUNK
 7189       FORMAT('      in the vicinity of ',F12.3)
            WRITE(IPR,999)
          ENDIF
        ENDIF
C
        ATEMP=REAL(NTOTZ)*THRPER + 0.00001
        MINTHR=INT(ATEMP)
        IF(NABOVE.GE.MINTHR)THEN
          IF(NFINAL.EQ.0)THEN
            NFINAL=I
          ELSE
            NFINAL=MIN(I,NFINAL)
          ENDIF
        ELSE
          NFINAL=0
        ENDIF
        IF(NN(I).GT.0)THEN
          YM(I)=REAL(DSUM/DBLE(NN(I)))
          CALL SORT(TEMP1,NPLUS,TEMP1,IPR)
        ELSE
          YM(I)=0.0
        ENDIF
        IF(I.EQ.1)THEN
          IF(YMINMS.LE.0.0)THEN
            HH=YMIN
          ELSE
            HH=2.0*YMIN - YMIN2
          ENDIF
C
          NUMPAD=NTOTZ-NPLUS
          ICNT=0
          IF(NUMPAD.GT.0)THEN
            DO1112II=1,NUMPAD
              ICNT=ICNT+1
              TEMP3(ICNT)=0.0
 1112       CONTINUE
          ENDIF
          DO1113II=1,NPLUS
            ICNT=ICNT+1
            TEMP3(ICNT)=TEMP1(II)
 1113     CONTINUE
C
CCCCC     P100=95.0
          P100=(1.0 - ALPHA)*100.
          CALL PERCEN(P100,TEMP3,NTOTZ,IPR,IQUAME,TEMP2,MAXOBV,CV1,
     1                IDEBUG,IERROR)
          IF(IERROR.EQ.'YES')GOTO9000
        ENDIF
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='3.1'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,1161)
 1161     FORMAT('***** DPLOD2: AFTER COMPUTE SAMPLE SIZES')
          WRITE(IPR,1162)I,MU(I),NN(I),NS(I),YM(I)
 1162     FORMAT('I,MU(I),NN(I),NS(I),YM(I) = ',4I8,G15.7)
          WRITE(IPR,1163)HH,THRPER,NABOVE,NFINAL,MINTHR,CV
 1163     FORMAT('HH,THRPER,NABOVE,NFINAL,MINTHR,CV = ',
     1           2F10.5,3I8,G15.7)
        ENDIF
C
C       STATISTICS FOR MASS = 0 VALUES
C
        IF(I.EQ.1)THEN
C
C         IF TOO FEW POSITIVE VALUES, THEN CLEAN CASE IS CALLED.
C
          YM1=YM(I)
          CALL VAR(TEMP1,NPLUS,IPR,VYN,IDEBUG,IERROR)
          NU1=NN(I)-1
          YV=REAL(NU1)*VYN/REAL(NN(I))
          CALL NORPPF((REAL(MU(I))+0.5)/(REAL(NS(I))+1.0),ZZ,IPR)
          QQ=MAX(ZA,ZZ)
          FT=0.5*ZZ*(YM(1)-HH)
          SIG0=FT + SQRT((FT**2) + YV + (YM(I)-HH)**2)
          CV2=YM1 + QQ*SIG0
          AF=HH - SIG0*ZZ
          IF(NU1.GE.1)THEN
            CALL CHSPPF(GA,NU1,GB1,IPR)
            GB=SQRT(REAL(NU1)/GB1)
            SU1=GB*SIG0
          ELSE
            GB=CPUMIN
            SU1=CPUMIN
          ENDIF
          SN=((QQ*SIG0)**2)/REAL(NN(I))
          IF(CVUSER.GT.0.0)THEN
            CV=CVUSER
          ELSE
            IF(ILODCV.EQ.'QUAN')THEN
              CV=CV1
            ELSE
              CV=CV2
            ENDIF
          ENDIF
C
          IF(IDEBUG.EQ.'ON')THEN
            ISTEPN='3.1'
            WRITE(IPR,991)ISTEPN
            WRITE(IPR,999)
            WRITE(IPR,1171)
 1171       FORMAT('***** DPLOD2: MASS = 0 CASE')
            WRITE(IPR,1172)HH,YM(1),VYN,YV,ZZ
 1172       FORMAT('HH,YM(1),VYN,YV,ZZ = ',5G15.7)
            WRITE(IPR,1173)QQ,FT,SIG0,AF,GA
 1173       FORMAT('QQ,FT,SIG0,AF,GA = ',5G15.7)
            WRITE(IPR,1174)NU1,GB1,GB,SU1,SN
 1174       FORMAT('NU1,GB1,GB,SU1,SN = ',I8,4G15.7)
          ENDIF
C
        ENDIF
C
C       COMPUTE MEAN OF Y - AF FOR NON-ZERO VALUES
C
        IF(NPLUS.GE.1)THEN
          CALL SD(TEMP1,NPLUS,IPR,XVAR,IDEBUG,IERROR)
        ELSE
          XVAR=0.0
        ENDIF
        YVAR(I)=XVAR
        SUM=0.0
        DO1180J=1,NPLUS
          TEMP1(J)=TEMP1(J) - AF
          SUM=SUM + TEMP1(J)**2
 1180   CONTINUE
        IF(NPLUS.GE.1)THEN
          CALL MEAN(TEMP1,NPLUS,IPR,YTTM,IDEBUG,IERROR)
        ELSE
          YTTM=0.0
        ENDIF
        YT(I)=YTTM
        YTT(I)=SUM
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='3.3'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,1192)I,YT(I),YTT(I)
 1192     FORMAT('I,YT(I),YTT(I) = ',I8,2G15.7)
        ENDIF
C
 1100 CONTINUE
C
      DO2000J=2,NDIST
C
        XXNN=0.0
        BAK=0.0
        SYY2=0.0
        ITN=0
        DO2010I=2,J
          XXNN=XXNN + XIDTEM(I)**2*REAL(NN(I))
          BAK=BAK + XIDTEM(I)*YT(I)*REAL(NN(I))
          SYY2=SYY2 + YTT(I)
          ITN=ITN + NN(I)
 2010   CONTINUE
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='4.1'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,2061)
 2061     FORMAT('***** DPLOD2: AT 2010')
          WRITE(IPR,2062)J,XXNN,BAK,SYY2,ITN
 2062     FORMAT('J,XXNN,BAK,SYY2,ITN = ',I8,3G15.7,I8)
          WRITE(IPR,2063)J,SIG0,ZB,QQ,GA
 2063     FORMAT('J,SIG0,ZB,QQ,GA = ',I8,4G15.7)
        ENDIF
C
CCCCC   MAY 2009: IF B1K1 IS NEGATIVE, THEN VALUES COMPUTED
CCCCC             FOR LU(J) AND TO(J) TEND TO BE NONSENSE.
CCCCC             SO IN THIS CASE, SET THEM TO "CPUMIN" (THESE
CCCCC             WILL BE PRINTED AS "**").
CCCCC
CCCCC   JULY 2009: IF GB (COMPUTED ABOVE) = CPUMIN, THEN LIKEWISE
CCCCC              SET LU(J) AND TO(J) TO CPUMIN.
C
        ITN=ITN-1
        IF(GB.EQ.CPUMIN .OR. ITN.LT.1)THEN
          BA(J)=CPUMIN
          SI(J)=CPUMIN
          LO(J)=CPUMIN
          UN(J)=CPUMIN
          LU(J)=CPUMIN
          TO(J)=CPUMIN
        ELSE
          DD=1.0/SQRT(XXNN)
          BA(J)=BAK/XXNN
          UA(J)=(SYY2 - BA(J)*BAK)/REAL(ITN)
          SI(J)=SQRT(UA(J))
          CALL TPPF(1.0 - 0.5*GA,REAL(ITN),TTK,IPR)
          T(J)=TTK
          B1K1=BA(J)-(DD*T(J)*SI(J))
          B1K2=0.000001
          B1(J)=MAX(B1K1,0.000001)
          CALL CHSPPF(GA/2.0,ITN,PPK,IPR)
          P(J)=PPK
          SU(J)=SI(J)*SQRT(REAL(ITN)/P(J))
          LO(J)=(QQ*SIG0 + ZB*SI(J))/BA(J)
          LU(J)=(QQ*SU1 + ZB*SU(J))/B1(J)
          S(J)=(SN + (AZ*UA(J)/REAL(ITN)) + UA(J)*(DD**2)*
     1         (LO(J)**2))/BA(J)/BA(J)
          UN(J)=SQRT(S(J))
          TO(J)=(QQ*SU1 + (ZB+ZP)*SU(J))/B1(J)
C
          IF(B1(J).EQ.B1K2)THEN
            LU(J)=CPUMIN
            TO(J)=CPUMIN
          ENDIF
C
        ENDIF
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='4.2'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,2072)J,DD,BA(J),UA(J),SI(J),TTK,B1K1
 2072     FORMAT('J,DD,BA(J),UA(J),SI(J),TTK,B1K1 = ',I8,6G15.7)
          WRITE(IPR,2073)J,B1(J),P(J),SU(J),LO(J)
 2073     FORMAT('J,B1(J),P(J),SU(J),LO(J) = ',I8,4G15.7)
          WRITE(IPR,2074)J,LU(J),S(J),UN(J),TO(J)
 2074     FORMAT('J,LU(J),S(J),UN(J),TO(J) = ',I8,4G15.7)
        ENDIF
C
 2000 CONTINUE
C
      SI(1)=SIG0
C
C     NOTE 4/2009: CHECK IF MORE BLANK MEASUREMENTS NEEDED.
C
C     NOTE 5/2009: THIS RULE CHANGED SO THAT A MAXIMUM OF 20
C                  BLANKS CALLED FOR ABOVE.
C
CCCCC DELTA=0.8
CCCCC M0=NN(1)
CCCCC CALL CHSPPF(DELTA,M0,ATEMP)
CCCCC STEMP=(SIG0/ATEMP) - 1.0
CCCCC IF(STEMP.LE.0.0)THEN
CCCCC   IFLAG9=0
CCCCC   NFINAL=NFINAL+1
CCCCC   IF(NFINAL.GT.NDIST)NFINAL=NDIST
CCCCC   ALOD=LO(NFINAL)
CCCCC   ALODSE=UN(NFINAL)
CCCCC ELSE
CCCCC   IFLAG9=1
CCCCC   AMORE=REAL(M0)*STEMP
CCCCC   NMORE=INT(AMORE+0.99)
CCCCC ENDIF
C
      IF(NFINAL.EQ.0)NFINAL=NDIST
      IFLAG9=0
      ALOD=LO(NFINAL)
      ALODSE=UN(NFINAL)
      ALODLU=LU(NFINAL)
C
      IF(IDEBUG.EQ.'ON')THEN
        ISTEPN='4.3'
        WRITE(IPR,991)ISTEPN
        WRITE(IPR,999)
        WRITE(IPR,2082)NFINAL
 2082   FORMAT('NFINAL,ALOD,ALODSE = ',I8,2G15.7)
      ENDIF
C
      IVAL1=INT(100.0*(1.0 - GAMMA)+0.1)
      IVAL2=INT(100.0*(1.0 - GAMMA)+0.1)
      IVAL3=INT(100.0*(1.0 - PA)+0.1)
C
C               **************************************************
C               **  STEP 5--                                    **
C               **   OPEN THE STORAGE FILES, WRITE RELEVANT     **
C               **   QUANTITIES TO FILE, AND THEN CLOSE         **
C               **   STORAGE FILES.                             **
C               **************************************************
C
      ISTEPN='5'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
C
      IF(ILODOF.EQ.'ON')THEN
        IOUNI1=11
        IOUNI2=12
        OPEN(IOUNI1,FILE='TABLE1.OUT')
        OPEN(IOUNI2,FILE='TABLE2.OUT')
        DO3010I=1,NDIST
          WRITE(IOUNI1,3012)XIDTEM(I),REAL(MU(I)),REAL(NN(I)),YT(I),
     1                      YVAR(I)
          WRITE(IOUNI2,3012)XIDTEM(I),BA(I),SI(I),LO(I),UN(I),LU(I),
     1                      TO(I)
 3012     FORMAT(7E15.7)
 3010   CONTINUE
        CLOSE(IOUNI1)
        CLOSE(IOUNI2)
      ENDIF
C
      IF(YCONST.GT.0.5)THEN
        CV=CV - YCONST
      ENDIF
C
C     WRITE "FINAL" ESTIMATES SEPARATELY
C
      IXX=NUMDIG
      IYY=NUMDIG+8
      IFORMT='(A37,Gyy.x)'
      WRITE(IFORMT(7:8),'(I2)')IYY
      WRITE(IFORMT(10:10),'(I1)')IXX
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,*)'IFORMT = ',IFORMT
      ENDIF
C
      WRITE(IPR,999)
      WRITE(IPR,999)
      WRITE(IPR,7090)
 7090 FORMAT('            LIMITS OF DETECTION ANALYSIS')
      WRITE(IPR,999)
      IF(IFLAG9.EQ.0)THEN
        WRITE(IPR,7091)
 7091   FORMAT('Final Estimate:')
        IF(ILODPC.EQ.'ON')THEN
          IF(ALPHA.EQ.0.05)THEN
            NCHAR(1)=37
            IHEAD(1:NCHAR(1))='Critical Value (CV95)              = '
          ELSEIF(ALPHA.EQ.0.10)THEN
             NCHAR(1)=37
             IHEAD(1:NCHAR(1))='Critical Value (CV90)              = '
          ELSEIF(ALPHA.EQ.0.01)THEN
             NCHAR(1)=37
             IHEAD(1:NCHAR(1))='Critical Value (CV99)              = '
          ELSE
            NCHAR(1)=37
            IHEAD(1:NCHAR(1))='Critical Value (CVxx.x)            = '
            ATEMP=100.*(1.0 - ALPHA)
            WRITE(IHEAD(19:22),'(F4.1)')ATEMP
          ENDIF
C
          IF(IDEBUG.EQ.'ON')THEN
            WRITE(IPR,*)'CV,NCHAR(1),IHEAD = ',CV,NCHAR(1),IHEAD
          ENDIF
C
          WRITE(IPR,IFORMT)IHEAD(1:NCHAR(1)),CV
        ENDIF
C
        IF(ALPHA.EQ.0.05)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Detection Limit (LOD95)            = '
        ELSEIF(ALPHA.EQ.0.10)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Detection Limit (LOD90)            = '
        ELSEIF(ALPHA.EQ.0.01)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Detection Limit (LOD99)            = '
        ELSE
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Detection Limit (LODxx.x)          = '
          ATEMP=100.*(1.0 - ALPHA)
          WRITE(IHEAD(20:23),'(F4.1)')ATEMP
        ENDIF
        WRITE(IPR,IFORMT)IHEAD(1:NCHAR(1)),ALOD
C
        IF(GAMMA.EQ.0.05)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='95% Upper Confidence Limit on LOD  = '
        ELSEIF(GAMMA.EQ.0.10)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='90% Upper Confidence Limit on LOD  = '
        ELSEIF(GAMMA.EQ.0.01)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='99% Upper Confidence Limit on LOD  = '
        ELSE
          IXX=NUMDIG
          IYY=NUMDIG+8
          IFORMT='(A38,Gyy.xx)'
          WRITE(IFORMT(7:8),'(I2)')IYY
          WRITE(IFORMT(10:11),'(I2)')IXX
          NCHAR(1)=38
          IHEAD(1:NCHAR(1))='xx.x% Upper Confidence Limit on LOD = '
          ATEMP=100.*(1.0 - GAMMA)
          WRITE(IHEAD(1:4),'(F4.1)')ATEMP
        ENDIF
        WRITE(IPR,IFORMT)IHEAD(1:NCHAR(1)),ALODLU
      ELSEIF(IFLAGZ)THEN
        WRITE(IPR,7081)
 7081   FORMAT('Measurements at each mass level are either')
        WRITE(IPR,7083)
 7083   FORMAT('all zero or all non-zero.  It is recommended')
        WRITE(IPR,7085)
 7085   FORMAT('that you take additional measurements in the')
        ATEMP1=(XIDTEM(IALLZ)+XIDTEM(IALLP))/2.0
        WRITE(IPR,7087)ATEMP1
 7087   FORMAT('vicinity of ',F12.3)
      ELSE
        WRITE(IPR,7093)NMORE
 7093   FORMAT('It is recommended that you take ',I5)
        WRITE(IPR,7095)
 7095   FORMAT('additional blank measurements.')
      ENDIF
      WRITE(IPR,999)
      WRITE(IPR,999)
      IHEAD=' '
      IFORMT=' '
      NCHAR(1)=0
C
C     PRINT TABLES
C
      IF(ILODST.EQ.'OFF')GOTO7199
C
      WRITE(IPR,7110)
 7110 FORMAT('Summary Table')
      WRITE(IPR,7112)
 7112 FORMAT('--------------------------------------------------',
     1       '------------------------------------')
C
      WRITE(IPR,7114)
 7114 FORMAT(11X,'Mass',6X,'Number of',10X,'Number of',12X,'Mean of',
     1       13X,'SD of')
      WRITE(IPR,7116)
 7116 FORMAT(9X,'Values',4X,'Zero Values',4X,'Non-Zero Values',
     1       4X,'Non-Zero Values',3X,'Non-Zero Values')
      WRITE(IPR,7112)
C
      DO7120I=1,NDIST
        WRITE(IPR,7121)XIDTEM(I),MU(I),NN(I),YM(I)-YCONST,YVAR(I)
 7121   FORMAT(F15.5,I15,I19,F19.5,F18.5)
 7120 CONTINUE
      WRITE(IPR,7112)
C
 7199 CONTINUE
C
      IF(ILODTA.EQ.'OFF')GOTO7299
C
      WRITE(IPR,7210)
 7210 FORMAT('LOD Table')
      WRITE(IPR,7212)
 7212 FORMAT('--------------------------------------------------',
     1       '-------------------------------------------------------')
C
      IVALC=' '
      IF(GAMMA.EQ.0.05.AND.PA.EQ.0.05)THEN
        IVALC='    95%-95%'
      ELSEIF(GAMMA.EQ.0.10.AND.PA.EQ.0.05)THEN
        IVALC='    90%-95%'
      ELSEIF(GAMMA.EQ.0.01.AND.PA.EQ.0.05)THEN
        IVALC='    99%-95%'
      ELSEIF(GAMMA.EQ.0.10.AND.PA.EQ.0.10)THEN
        IVALC='    90%-90%'
      ELSEIF(GAMMA.EQ.0.05.AND.PA.EQ.0.10)THEN
        IVALC='    95%-90%'
      ELSEIF(GAMMA.EQ.0.01.AND.PA.EQ.0.10)THEN
        IVALC='    99%-90%'
      ELSEIF(GAMMA.EQ.0.10.AND.PA.EQ.0.01)THEN
        IVALC='    90%-99%'
      ELSEIF(GAMMA.EQ.0.05.AND.PA.EQ.0.01)THEN
        IVALC='    95%-99%'
      ELSEIF(GAMMA.EQ.0.01.AND.PA.EQ.0.01)THEN
        IVALC='    99%-99%'
      ELSE
        IVALC='xx.x%-xx.x%'
        ATEMP=100.*(1.0 - GAMMA)
        WRITE(IVALC(1:4),'(F4.1)')ATEMP
        ATEMP=100.*(1.0 - PA)
        WRITE(IVALC(7:10),'(F4.1)')ATEMP
      ENDIF
      WRITE(IPR,7214)IVALC(1:11)
 7214 FORMAT(24X,'Linear',6X,'StdDev of',3X,'LOD Estimate',
     1       8X,'1-Sigma',19X,A11)
      IVALC=' '
      IF(GAMMA.EQ.0.05)THEN
        IVALC='      95'
      ELSEIF(GAMMA.EQ.0.10)THEN
        IVALC='      90'
      ELSEIF(GAMMA.EQ.0.01)THEN
        IVALC='      99'
      ELSE
        ATEMP=100.*(1.0 - GAMMA)
        IVALC='    xx.x'
        WRITE(IVALC(5:8),'(F4.1)')ATEMP
      ENDIF
      WRITE(IPR,7216)IVALC(1:8)
 7216 FORMAT(17X,'Least Squares',4X,'Errors From',7X,'Based On',
     1       4X,'Uncertainty',A8,'% Upper',10X,'Upper')
      WRITE(IPR,7218)
 7218 FORMAT(11X,'Mass',10X,'Slope',9X,'Linear',6X,'Cumulated',
     1       9X,'on LOD',5X,'Confidence',6X,'Tolerance')
      WRITE(IPR,7219)
 7219 FORMAT(9X,'Values',6X,'Estimates',10X,'Model',9X,'Sample',
     1       7X,'Estimate',3X,'Limit on LOD',3X,'Limit on LOD')
      WRITE(IPR,7212)
C
      WRITE(IPR,7231)XIDTEM(I)
 7231 FORMAT(F15.5,6(13X,'**'))
C
      DO7230I=2,NDIST
        WRITE(IPR,7232)XIDTEM(I),BA(I),SI(I),LO(I),UN(I),LU(I),TO(I)
 7232   FORMAT(7F15.5)
 7230 CONTINUE
C
 7299 CONTINUE
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT       **
C               *****************
C
 9000 CONTINUE
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF DPLOD2--')
        WRITE(IPR,9012)IERROR
 9012   FORMAT('IERROR = ',A4)
      ENDIF
C
      RETURN
      END
      SUBROUTINE DPLOD3(Y,X,N,MAXOBV,
     1                  XIDTEM,TEMP1,TEMP2,TEMP3,YT,YM,
     1                  D,BA,BL,SU,T,P,LO,LU,S,UN,TO,UA,SI,B1,YTT,YVAR,
     1                  NS,MU,NN,
     1                  ALPHA,BETA,GAMMA,PA,THRPER,CVUSER,
     1                  ALOD,ALODSE,CV,ALODLU,
     1                  YMINMS,YMIN,YMIN2,
     1                  ILODCV,ILODOF,ILODST,ILODTA,ILODPC,
     1                  IDEBUG,IERROR)
C
C     PURPOSE--PERFORM A LIMITS OF DETECTION ANALYIS.  SPECIFICALLY,
C              WE ARE CHECKING FOR A LOWER LIMIT OF DETECTION.
C              IT IS ASSUMED THAT WE CANNOT MEASURE THE MASS BELOW
C              A CERTAIN DETECTION LIMIT AND WE NEED TO DETERMINE
C              WHAT THAT DETECTION LIMIT IS.  THERE ARE 2 CASES:
C
C                  1) FOR A ZERO MASS, THE MEASUREMENT WILL ALWAYS
C                     RETURN A ZERO VALUE.  THIS WILL BE REFERRED
C                     TO AS THE "CLEAN" CASE".
C
C                  2) FOR A ZERO MASS, THE MEASUREMENT MAY BE A
C                     NON-ZERO VALUE.  THIS WILL BE REFERRED TO AS
C                     THE DIRTY CASE.
C
C              THE DATA CONSISTS OF A MASS VALUE ("X") AND A
C              MEASURED VALUE ("Y").  THE COMMAND SYNTAX IS:
C
C                  CLEAN LIMIT OF DETECTION  Y   X
C                  DIRTY LIMIT OF DETECTION  Y   X
C                  LIMIT OF DETECTION  Y   X
C
C              IF THE CLEAN/DIRTY CASE IS NOT SPECIFIED, ASSUME
C              THE DIRTY CASE.  THIS SUBROUTINE PERFORMS THE
C              COMPUTATIONS FOR THE "CLEAN" CASE.  PERFORM A CHECK
C              THAT ALL MEASUREMENTS FOR THE ZERO MASS ARE INDEED
C              ZERO (IF NOT, PRINT AN ERROR AND RETURN).
C
C     PRINTING--YES
C     SUBROUTINES NEEDED--NORPPF, CHSPPF
C     WRITTEN BY--ALAN HECKERT
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2899
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--2018/08
C     ORIGINAL VERSION--AUGUST    2018.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES--------------
C
      CHARACTER*4 ILODCV
      CHARACTER*4 ILODOF
      CHARACTER*4 ILODST
      CHARACTER*4 ILODTA
      CHARACTER*4 ILODPC
      CHARACTER*4 ISUBN0
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
      CHARACTER*4 IPTEMP
      CHARACTER*4 ISTEPN
C
      CHARACTER*40 IFORMT
      CHARACTER*60 IFRMT5
      CHARACTER*60 IFRMT6
      CHARACTER*15 ILABEL
      CHARACTER*11 IVALC
      CHARACTER*17 ITYPE
      CHARACTER*132 IHEAD
C
C----------------------------------------------------------------
C
      REAL Y(*)
      REAL X(*)
      REAL TEMP1(*)
      REAL TEMP2(*)
      REAL TEMP3(*)
      REAL YT(*)
      REAL YM(*)
      REAL XIDTEM(*)
C
      INTEGER NS(*)
      INTEGER MU(*)
      INTEGER NN(*)
C
      REAL D(*)
      REAL BA(*)
      REAL BL(*)
      REAL SU(*)
      REAL T(*)
      REAL P(*)
      REAL LO(*)
      REAL LU(*)
      REAL S(*)
      REAL UN(*)
      REAL TO(*)
      REAL UA(*)
      REAL SI(*)
      REAL B1(*)
      REAL YTT(*)
      REAL YVAR(*)
C
      REAL ZERO
      REAL EPS
C
      DOUBLE PRECISION DSUM
C
      INTEGER      NCHAR(8)
C
C
      LOGICAL IFLAGY
      LOGICAL IFLAGZ
C
      DATA ZERO /0.0/
      DATA EPS  /0.1E-7/
C
C-----START POINT------------------------------------------------
C
      IERROR='NO'
C
      NUMDIG=5
      IF(THRPER.EQ.CPUMIN)THRPER=0.6
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF DPLOD3--')
        WRITE(IPR,52)N,ALPHA,BETA,GAMMA,PA
   52   FORMAT('N,ALPHA,BETA,GAMMA,PA = ',I8,4G15.7)
        WRITE(IPR,53)THRPER,CVUSER
   53   FORMAT('THRPER,CVUSER = ',2G15.7)
        DO55I=1,N
          WRITE(IPR,56)I,Y(I),X(I)
   56     FORMAT('I,Y(I),X(I) = ',I8,2G15.7)
   55   CONTINUE
      ENDIF
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               **  BASICALLY, NEGATIVE VALUES ARE        **
C               **  CONSIDERED ERRORS.                    **
C               ********************************************
C
      IF(N.LT.2)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
  101   FORMAT('***** ERROR IN LIMITS OF DETECTION--')
        WRITE(IPR,102)
  102   FORMAT('      THE NUMBER OF OBSERVATIONS FOR THE LIMITS OF ',
     1         'DETECTION')
        WRITE(IPR,103)N
  103   FORMAT('      MUST BE AT LEAST 2;  THE ENTERED NUMBER OF ',
     1         'OBSERVATIONS = ',I8)
        WRITE(IPR,999)
        IERROR='YES'
        GOTO9000
      ENDIF
C
      DO200I=1,N
        IF(ABS(X(I)).LE.EPS)X(I)=ZERO
        IF(Y(I).LT.ZERO)THEN
          WRITE(IPR,999)
          WRITE(IPR,101)
          WRITE(IPR,202)
  202     FORMAT('      A NEGATIVE VALUE WAS ENCOUNTERED FOR ',
     1           'THE MEASURED VALUE')
          WRITE(IPR,204)I,Y(I)
  204     FORMAT('      ROW ',I8,' HAS THE VALUE ',G15.7)
          IERROR='YES'
          GOTO9000
        ELSEIF(X(I).LT.ZERO)THEN
          WRITE(IPR,999)
          WRITE(IPR,101)
          WRITE(IPR,212)
  212     FORMAT('      A NEGATIVE VALUE WAS ENCOUNTERED FOR ',
     1           'THE MASS VALUE')
          WRITE(IPR,204)I,X(I)
          IERROR='YES'
          GOTO9000
        ELSEIF(X(I).EQ.ZERO .AND. Y(I).GT.EPS)THEN
          WRITE(IPR,999)
          WRITE(IPR,101)
          WRITE(IPR,222)
  222     FORMAT('      FOR A ZERO MASS VALUE, A NON-ZERO ',
     1           'MEASUREMENT WAS DETECTED.')
          WRITE(IPR,223)
  223     FORMAT('      FOR THIS CASE, YOU NEED TO USE THE DIRTY ',
     1           'LIMITS OF DETECTION COMMAND')
          WRITE(IPR,224)
  224     FORMAT('      RATHER THAN THE CLEAN LIMITS OF DETECTION ',
     1           'COMMAND.')
          WRITE(IPR,228)I,X(I),Y(I)
  228     FORMAT('      ROW ',I8,' HAS MASS VALUE ',G15.7, ' AND ',
     1           'MEASUREMENT ',G15.7)
          IERROR='YES'
          GOTO9000
        ENDIF
  200 CONTINUE
C
C               ***********************************************
C               **  STEP 2.1--                               **
C               **  PERFORM THE BASIC CALCULATIONS.          **
C               ***********************************************
C
      ISTEPN='2.1'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
  991 FORMAT('AT STEP ',A4)
C
      CALL SORTC(X,Y,N,X,TEMP1,IPR)
      DO1005I=1,N
        Y(I)=TEMP1(I)
 1005 CONTINUE
C
C     ALPHA   = CRITICAL LEVEL
C     BETA    = LOD UNCERTAINTY CONFIDENCE
C     GAMMA   = UPPER CONFIDENCE BOUND
C     PA      = TOLERANCE BOUND SIGNIFICANCE
C
CCCCC NOTE 6/2010: ACCORDING TO ANDREW RUKHIN, THE ZA TERM SHOULD NOT BE
CCCCC              SQUARED.
C
      CALL NORPPF(1.0-ALPHA,ZA,IPR)
      CALL NORPPF(1.0-BETA,ZB,IPR)
      CALL NORPPF(1.0-PA,ZP,IPR)
      GA=GAMMA
C
      IF(IDEBUG.EQ.'ON')THEN
        ISTEPN='2.2'
        WRITE(IPR,991)ISTEPN
        WRITE(IPR,999)
        WRITE(IPR,61)
   61   FORMAT('***** FROM DPLOD3--')
        WRITE(IPR,62)ZA,AZ,ZB,ZP
   62   FORMAT('ZA,AZ,ZB,ZP = ',4G15.7)
      ENDIF
C
C     COMPUTE:
C
C       1)  NUMBER OF BLANK (ZERO) MASS VALUES           (M)
C       2)  NUMBER OF DISTINCT MASS VALUES               (NDIST)
C       3)  SAMPLE SIZES FOR EACH GROUP                  (NS)
C       4)  NUMBER OF ZERO RESPONSES IN EACH GROUP       (MU)
C       5)  NUMBER OF NON-ZERO RESPONSES IN EACH GROUP   (NN)
C       6)  MEAN OF NON-ZERO RESPONSES IN EACH GROUP     (YM)
C
C     FOR THE MASS = 0 CASE, COMPUTE SOME STATISTICS
C
      CALL DISTIN(X,N,XIDTEM,NDIST,IPR,IDEBUG,IERROR)
C
      IF(NDIST.LT.2)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
        WRITE(IPR,1012)
 1012   FORMAT('      THE NUMBER OF DISTINCT MASS VALUES IS LESS ',
     1         'THAN 2')
        IERROR='YES'
        GOTO9000
      ELSEIF(NDIST.EQ.N)THEN
        WRITE(IPR,999)
        WRITE(IPR,101)
        WRITE(IPR,1013)
 1013   FORMAT('      THE NUMBER OF DISTINCT MASS VALUES IS EQUAL TO ',
     1         'THE SAMPLE SIZE')
        IERROR='YES'
        GOTO9000
      ENDIF
C
      NFINAL=0
      IFLAGY=.FALSE.
      IFLAGZ=.TRUE.
      IALLZ=-99
      IALLP=-99
      ISTRT=2
C
      DO1100I=1,NDIST
C
C       NOTE: IF THERE ARE "ZERO" MASS VALUES (I.E., BLANKS),
C             THEN THESE NEED TO BE HANDLED SEPARATELY.
C
        HOLD=XIDTEM(I)
C
        D(I)=0.0
        BL(I)=0.0
        SU(I)=0.0
        T(I)=0.0
        P(I)=0.0
        IF(I.GT.1 .OR. (I.EQ.1 .AND. HOLD.EQ.ZERO))THEN
          LO(I)=0.0
          LU(I)=0.0
          BA(I)=0.0
          UN(I)=0.0
          S(I)=0.0
          TO(I)=0.0
        ELSE
          LO(I)=99.0
          LU(I)=99.0
          BA(I)=99.0
          UN(I)=99.0
          S(I)=99.0
          TO(I)=99.0
        ENDIF
C
        NZERO=0
        NPLUS=0
        DSUM=0.0D0
        DO1110J=1,N
          IF(X(J).EQ.HOLD)THEN
            IF(Y(J).LE.EPS)THEN
              NZERO=NZERO+1
            ELSE
              NPLUS=NPLUS+1
              TEMP1(NPLUS)=Y(J)
              DSUM=DSUM + DBLE(Y(J))
            ENDIF
          ENDIF
 1110   CONTINUE
        NTOTZ=NZERO+NPLUS
        ATEMP=REAL(NTOTZ)*(1.0 - THRPER)+0.00001
        MINZER=INT(ATEMP)
        MU(I)=NZERO
        NN(I)=NPLUS
        NS(I)=NTOTZ
C
        IF(NZERO.GE.1 .AND. NPLUS.GE.1)IFLAGZ=.FALSE.
        IF(IFLAGZ)THEN
          IF(NZERO.EQ.NTOTZ)IALLZ=I
          IF(NPLUS.EQ.NTOTZ .AND. IALLP.LT.0)IALLP=I
        ENDIF
        IF(I.EQ.NDIST)THEN
          NTEMPZ=INT(0.90*REAL(NTOTZ)+0.1)
          IF(NZERO.GT.NTEMPZ)THEN
            IFLAGY=.TRUE.
            WRITE(IPR,999)
            WRITE(IPR,101)
            WRITE(IPR,7071)
 7071       FORMAT('      More than 90% of the measurements at the')
            WRITE(IPR,7073)
 7073       FORMAT('      highest mass level are zero.')
            WRITE(IPR,7075)
 7075       FORMAT('      The LOD analysis will not be performed.')
            WRITE(IPR,999)
            IERROR='YES'
            GOTO9000
          ENDIF
        ENDIF
C
        IF(NZERO.GT.MINZER)NFINAL=I
        IF(NN(I).GT.0)THEN
          YM(I)=REAL(DSUM/DBLE(NN(I)))
          CALL SORT(TEMP1,NPLUS,TEMP1,IPR)
        ELSE
          YM(I)=0.0
        ENDIF
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='3.1'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,1061)
 1061     FORMAT('***** DPLOD3: AFTER COMPUTE SAMPLE SIZES')
          WRITE(IPR,1062)I,MU(I),NN(I),NS(I),YM(I)
 1062     FORMAT('I,MU(I),NN(I),NS(I),YM(I) = ',4I8,G15.7)
          WRITE(IPR,1063)NZERO,MINZER,NFINAL
 1063     FORMAT('NZERO,MINZER,NFINAL = ',3I8)
        ENDIF
C
C       STATISTICS FOR MASS = 0 VALUES
C
C       COMPUTE SUM OF Y**2 FOR NON-ZERO VALUES
C
        IF(NPLUS.GE.1)THEN
          CALL SD(TEMP1,NPLUS,IPR,XVAR,IDEBUG,IERROR)
          CALL MEAN(TEMP1,NPLUS,IPR,XMEAN,IDEBUG,IERROR)
        ELSE
          XVAR=0.0
          XMEAN=0.0
        ENDIF
        YM(I)=XMEAN
        YVAR(I)=XVAR
        SUM=0.0
        DO1030J=1,NPLUS
          SUM=SUM + TEMP1(J)**2
 1030   CONTINUE
        YTT(I)=SUM
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='3.3'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,1082)I,YTT(I)
 1082     FORMAT('I,YTT(I) = ',I8,G15.7)
        ENDIF
C
 1100 CONTINUE
C
C
CCCCC   APRIL 2010: X1 SHOULD BE SET TO X(1) IF ALL
CCCCC               NON-ZERO MASSES HAVE AT LEAST ONE NON-ZERO
CCCCC               RESPONSE.  HOWEVER, IF SOME NON-ZERO MASSES
CCCCC               HAVE ALL ZERO RESPONSES, THEN SET X1 TO X(K)
CCCCC               WHERE K IS THE SMALLEST MASS LEVEL WITH AT LEAST
CCCCC               ONE NON-ZERO RESPONSE.
CCCCC
CCCCC   MAY 2010:   NOTE THAT DEGREES OF FREEDOM FOR CHI-SQUARE AND
CCCCC               T PERCENT POINT ROUTINES BELOW REQUIRES THAT CUMULATIVE
CCCCC               NUMBER OF NON-ZERO VALUES BE AT LEAST TWO.  SO
CCCCC               ACTUALLY LOOK FOR FIRST MASS LEVEL WITH AT LEAST TWO
CCCCC               NON-ZERO RESPONSES.
CCCCC
CCCCC   MAY 2010:   ACTUALLY NEED TO SET X1 TO X(K) WHERE K IS THE
CCCCC               THE SMALLEST MASS LEVEL WITH A NON-ZERO VARIANCE
CCCCC               (I.E., IF HAVE 2 POINTS BUT THEY HAVE THE SAME VALUE,
CCCCC               CANNOT COMPUTE A NON-ZERO VARIANCE).
CCCCC
CCCCC X1=X(1)
      KTEMP=1
      DO1907II=2,NDIST
        IF(YVAR(II).GT.0.0)THEN
          KTEMP=II
          GOTO1908
        ENDIF
 1907 CONTINUE
 1908 CONTINUE
      X1=XIDTEM(KTEMP)
      IF(KTEMP.GT.ISTRT)ISTRT=KTEMP
C
      CALL NORPPF((REAL(MU(KTEMP))+0.5)/(REAL(NS(KTEMP))+1.0),ZZ,IPR)
      QQ=MAX(ZA,ZZ)
C
C     COMPUTE CRITICAL VALUE FOR CLEAN CASE.
C
C     NOTE: CLARIFICATION FROM MICHAEL VERKOUTEREN - FOR CLEAN CASE,
C           THE REPORTED CRITICAL VALUE IS THE "THRESHOLD" (THE HH
C           TERM BELOW), NOT THE 95-TH PERCENTILE OF THE FIRST MASS
C           LEVEL WHERE WE HAVE DATA.
C
      KTEMP2=KTEMP
      IF(KTEMP.EQ.1)KTEMP2=2
      IF(CVUSER.GT.0.0)THEN
        HH=CVUSER
      ELSE
        HH=2.0*YMIN - YMIN2
      ENDIF
C
      DO2000J=ISTRT,NDIST
C
        XXNN=0.0
        BAK=0.0
        SYY2=0.0
        ITN=0
        DO2010I=ISTRT,J
          XXNN=XXNN + XIDTEM(I)**2*REAL(NN(I))
          BAK=BAK + XIDTEM(I)*YM(I)*REAL(NN(I))
          SYY2=SYY2 + YTT(I)
          ITN=ITN + NN(I)
 2010   CONTINUE
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='4.1'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,2061)
 2061     FORMAT('***** DPLOD3: AT 2010')
          WRITE(IPR,2062)J,XXNN,BAK,SYY2,ITN
 2062     FORMAT('J,XXNN,BAK,SYY2,ITN = ',I8,3G15.7,I8)
          WRITE(IPR,2063)ZA,ZB,QQ,X1
 2063     FORMAT('ZA,ZB,QQ,X1 = ',4G15.7)
        ENDIF
C
        DD=1.0/SQRT(XXNN)
        BA(J)=BAK/XXNN
        ITN=ITN-1
        UA(J)=(SYY2 - BA(J)*BAK)/REAL(ITN)
        SI(J)=SQRT(UA(J))
        CALL TPPF(1.0 - 0.5*GA,REAL(ITN),TTK)
        T(J)=TTK
        B1K1=BA(J)-(DD*T(J)*SI(J))
        B1K2=0.000001
        B1(J)=MAX(B1K1,0.000001)
        CALL CHSPPF(GA/2.0,ITN,PPK,IPR)
        P(J)=PPK
        SU(J)=SI(J)*SQRT(REAL(ITN)/P(J))
        LO(J)=X1 + ((QQ + ZB)*SI(J))/BA(J)
        LU(J)=X1 + ((QQ+ZB)*SU(J))/B1(J)
        S(J)=((AZ*UA(J)/REAL(ITN)) + UA(J)*(DD**2)*(LO(J)**2))/
     1       BA(J)/BA(J)
        UN(J)=SQRT(S(J))
        TO(J)=X1 + ((QQ+ZB+ZP)*SU(J))/B1(J)
C
        IF(IDEBUG.EQ.'ON')THEN
          ISTEPN='4.2'
          WRITE(IPR,991)ISTEPN
          WRITE(IPR,999)
          WRITE(IPR,2072)J,DD,BA(J),UA(J),SI(J),TTK
 2072     FORMAT('J,DD,BA(J),UA(J),SI(J),TTK = ',I8,5G15.7)
          WRITE(IPR,2073)J,B1(J),P(J),SU(J),LO(J)
 2073     FORMAT('J,B1(J),P(J),SU(J),LO(J) = ',I8,4G15.7)
          WRITE(IPR,2074)J,LU(J),S(J),UN(J),TO(J)
 2074     FORMAT('J,LU(J),S(J),UN(J),TO(J) = ',I8,4G15.7)
        ENDIF
C
 2000 CONTINUE
C
      SI(1)=SIG0
C
      NFINAL=NFINAL+1
      IF(NFINAL.GT.NDIST)NFINAL=NDIST
      ALOD=LO(NFINAL)
      ALODSE=UN(NFINAL)
      ALODLU=LU(NFINAL)
C
      IF(IDEBUG.EQ.'ON')THEN
        ISTEPN='4.3'
        WRITE(IPR,991)ISTEPN
        WRITE(IPR,999)
        WRITE(IPR,2082)NFINAL
 2082   FORMAT('NFINAL,ALOD,ALODSE = ',I8,2G15.7)
      ENDIF
C
      IVAL1=INT(100.0*(1.0 - GAMMA)+0.1)
      IVAL2=INT(100.0*(1.0 - GAMMA)+0.1)
      IVAL3=INT(100.0*(1.0 - PA)+0.1)
C
C               **************************************************
C               **  STEP 5--                                    **
C               **   OPEN THE STORAGE FILES, WRITE RELEVANT     **
C               **   QUANTITIES TO FILE, AND THEN CLOSE         **
C               **   STORAGE FILES.                             **
C               **************************************************
C
      ISTEPN='5'
      IF(IDEBUG.EQ.'ON')WRITE(IPR,991)ISTEPN
C
      IF(ILODOF.EQ.'ON')THEN
        IOUNI1=11
        IOUNI2=12
        OPEN(IOUNI1,FILE='TABLE1.OUT')
        OPEN(IOUNI2,FILE='TABLE2.OUT')
        DO3010I=1,NDIST
          WRITE(IOUNI1,3012)XIDTEM(I),REAL(MU(I)),REAL(NN(I)),YT(I),
     1                      YVAR(I)
          WRITE(IOUNI2,3012)XIDTEM(I),BA(I),SI(I),LO(I),UN(I),LU(I),
     1                      TO(I)
 3012     FORMAT(7E15.7)
 3010   CONTINUE
        CLOSE(IOUNI1)
        CLOSE(IOUNI2)
      ENDIF
C
      CV=HH
C
C     WRITE "FINAL" ESTIMATES SEPARATELY
C
C     2014/03: TO BE CONSISTENT WITH ASTM E2677 STANDARD, CHANGE
C              "Scxx.x"  to   "CVxx.x" AND CHANGE  "LDxx"  TO  "LODxx".
C
      IXX=NUMDIG
      IYY=NUMDIG+8
      IFORMT='(A37,Gyy.xx)'
      WRITE(IFORMT(7:8),'(I2)')IYY
      WRITE(IFORMT(10:11),'(I2)')IXX
C
      WRITE(IPR,999)
      WRITE(IPR,999)
      WRITE(IPR,7090)
 7090 FORMAT('            LIMITS OF DETECTION ANALYSIS')
      WRITE(IPR,999)
      IF(IFLAGZ)THEN
        WRITE(IPR,7081)
 7081   FORMAT('Measurements at each mass level are either')
        WRITE(IPR,7083)
 7083   FORMAT('all zero or all non-zero.  It is recommended')
        WRITE(IPR,7085)
 7085   FORMAT('that you take additional measurements in the')
        ATEMP1=(XIDTEM(IALLZ)+XIDTEM(IALLP))/2.0
        WRITE(IPR,7087)ATEMP1
 7087   FORMAT('vicinity of ',F12.3)
        WRITE(IPR,999)
        WRITE(IPR,999)
      ENDIF
C
      WRITE(IPR,7091)
 7091 FORMAT('Final Estimate:')
      IF(ILODPC.EQ.'ON')THEN
        IF(ALPHA.EQ.0.05)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Critical Value (CV95)              = '
        ELSEIF(ALPHA.EQ.0.10)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Critical Value (CV90)              = '
        ELSEIF(ALPHA.EQ.0.01)THEN
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Critical Value (CV99)              = '
        ELSE
          NCHAR(1)=37
          IHEAD(1:NCHAR(1))='Critical Value (CVxx.x)            = '
          ATEMP=100.*(1.0 - ALPHA)
          WRITE(IHEAD(19:22),'(F4.1)')ATEMP
        ENDIF
        WRITE(IPR,IFORMT)IHEAD(1:NCHAR(1)),HH
      ENDIF
C
      IF(ALPHA.EQ.0.05)THEN
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='Detection Limit (LOD95)            = '
      ELSEIF(ALPHA.EQ.0.10)THEN
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='Detection Limit (LOD90)            = '
      ELSEIF(ALPHA.EQ.0.01)THEN
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='Detection Limit (LOD99)            = '
      ELSE
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='Detection Limit (LODxx.x)          = '
        ATEMP=100.*(1.0 - ALPHA)
        WRITE(IHEAD(20:23),'(F4.1)')ATEMP
      ENDIF
      WRITE(IPR,IFORMT)IHEAD(1:NCHAR(1)),ALOD
      IF(GAMMA.EQ.0.05)THEN
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='95% Upper Confidence Limit on LOD  = '
      ELSEIF(GAMMA.EQ.0.10)THEN
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='90% Upper Confidence Limit on LOD  = '
      ELSEIF(GAMMA.EQ.0.01)THEN
        NCHAR(1)=37
        IHEAD(1:NCHAR(1))='99% Upper Confidence Limit on LOD  = '
      ELSE
        IXX=NUMDIG
        IYY=NUMDIG+8
        IFORMT='(A38,Gyy.xx)'
        WRITE(IFORMT(7:8),'(I2)')IYY
        WRITE(IFORMT(10:11),'(I2)')IXX
        NCHAR(1)=38
        IHEAD(1:NCHAR(1))='xx.x% Upper Confidence Limit on LOD = '
        ATEMP=100.*(1.0 - GAMMA)
        WRITE(IHEAD(1:4),'(F4.1)')ATEMP
      ENDIF
      WRITE(IPR,IFORMT)IHEAD(1:NCHAR(1)),ALODLU
C
      WRITE(IPR,999)
      WRITE(IPR,999)
      IHEAD=' '
      IFORMT=' '
      NCHAR(1)=0
C
C     PRINT TABLES
C
      IF(ILODST.EQ.'OFF')GOTO7199
C
      WRITE(IPR,7110)
 7110 FORMAT('Summary Table')
      WRITE(IPR,7112)
 7112 FORMAT('--------------------------------------------------',
     1       '------------------------------------')
C
      WRITE(IPR,7114)
 7114 FORMAT(11X,'Mass',6X,'Number of',10X,'Number of',12X,'Mean of',
     1       13X,'SD of')
      WRITE(IPR,7116)
 7116 FORMAT(9X,'Values',4X,'Zero Values',4X,'Non-Zero Values',
     1       4X,'Non-Zero Values',3X,'Non-Zero Values')
      WRITE(IPR,7112)
C
      DO7120I=1,NDIST
        WRITE(IPR,7121)XIDTEM(I),MU(I),NN(I),YM(I)-YCONST,YVAR(I)
 7121   FORMAT(F15.5,I15,I19,F19.5,F18.5)
 7120 CONTINUE
      WRITE(IPR,7112)
C
 7199 CONTINUE
C
      IF(ILODTA.EQ.'OFF')GOTO7299
C
      WRITE(IPR,7210)
 7210 FORMAT('LOD Table')
      WRITE(IPR,7212)
 7212 FORMAT('--------------------------------------------------',
     1       '-------------------------------------------------------')
C
      IVALC=' '
      IF(GAMMA.EQ.0.05.AND.PA.EQ.0.05)THEN
        IVALC='    95%-95%'
      ELSEIF(GAMMA.EQ.0.10.AND.PA.EQ.0.05)THEN
        IVALC='    90%-95%'
      ELSEIF(GAMMA.EQ.0.01.AND.PA.EQ.0.05)THEN
        IVALC='    99%-95%'
      ELSEIF(GAMMA.EQ.0.10.AND.PA.EQ.0.10)THEN
        IVALC='    90%-90%'
      ELSEIF(GAMMA.EQ.0.05.AND.PA.EQ.0.10)THEN
        IVALC='    95%-90%'
      ELSEIF(GAMMA.EQ.0.01.AND.PA.EQ.0.10)THEN
        IVALC='    99%-90%'
      ELSEIF(GAMMA.EQ.0.10.AND.PA.EQ.0.01)THEN
        IVALC='    90%-99%'
      ELSEIF(GAMMA.EQ.0.05.AND.PA.EQ.0.01)THEN
        IVALC='    95%-99%'
      ELSEIF(GAMMA.EQ.0.01.AND.PA.EQ.0.01)THEN
        IVALC='    99%-99%'
      ELSE
        IVALC='xx.x%-xx.x%'
        ATEMP=100.*(1.0 - GAMMA)
        WRITE(IVALC(1:4),'(F4.1)')ATEMP
        ATEMP=100.*(1.0 - PA)
        WRITE(IVALC(7:10),'(F4.1)')ATEMP
      ENDIF
      WRITE(IPR,7214)IVALC(1:11)
 7214 FORMAT(24X,'Linear',6X,'StdDev of',3X,'LOD Estimate',
     1       8X,'1-Sigma',19X,A11)
      IVALC=' '
      IF(GAMMA.EQ.0.05)THEN
        IVALC='      95'
      ELSEIF(GAMMA.EQ.0.10)THEN
        IVALC='      90'
      ELSEIF(GAMMA.EQ.0.01)THEN
        IVALC='      99'
      ELSE
        ATEMP=100.*(1.0 - GAMMA)
        IVALC='    xx.x'
        WRITE(IVALC(5:8),'(F4.1)')ATEMP
      ENDIF
      WRITE(IPR,7216)IVALC(1:8)
 7216 FORMAT(17X,'Least Squares',4X,'Errors From',7X,'Based On',
     1       4X,'Uncertainty',A8,'% Upper',10X,'Upper')
      WRITE(IPR,7218)
 7218 FORMAT(11X,'Mass',10X,'Slope',9X,'Linear',6X,'Cumulated',
     1       9X,'on LOD',5X,'Confidence',6X,'Tolerance')
      WRITE(IPR,7219)
 7219 FORMAT(9X,'Values',6X,'Estimates',10X,'Model',9X,'Sample',
     1       7X,'Estimate',3X,'Limit on LOD',3X,'Limit on LOD')
      WRITE(IPR,7212)
C
      WRITE(IPR,7231)XIDTEM(I)
 7231 FORMAT(F15.5,6(13X,'**'))
C
      DO7230I=2,NDIST
        WRITE(IPR,7232)XIDTEM(I),BA(I),SI(I),LO(I),UN(I),LU(I),TO(I)
 7232   FORMAT(7F15.5)
 7230 CONTINUE
C
 7299 CONTINUE
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT       **
C               *****************
C
 9000 CONTINUE
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF DPLOD3--')
        WRITE(IPR,9012)IERROR
 9012   FORMAT('IERROR = ',A4)
      ENDIF
C
      RETURN
      END
      SUBROUTINE CHSPPF(P,NU,PPF,IPR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE FOR THE CHI-SQUARED DISTRIBUTION
C              WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU.
C              THE CHI-SQUARED DISTRIBUTION USED
C              HEREIN IS DEFINED FOR ALL NON-NEGATIVE X,
C              AND ITS PROBABILITY DENSITY FUNCTION IS GIVEN
C              IN REFERENCES 2, 3, AND 4 BELOW.
C              NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION
C              IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE
C              DISTRIBUTION FUNCTION OF THE DISTRIBUTION.
C     INPUT  ARGUMENTS--P      = THE SINGLE PRECISION VALUE
C                                (BETWEEN 0.0 (INCLUSIVELY)
C                                AND 1.0 (EXCLUSIVELY))
C                                AT WHICH THE PERCENT POINT
C                                FUNCTION IS TO BE EVALUATED.
C                     --NU     = THE INTEGER NUMBER OF DEGREES
C                                OF FREEDOM.
C                                NU SHOULD BE POSITIVE.
C     OUTPUT ARGUMENTS--PPF    = THE SINGLE PRECISION PERCENT
C                                POINT FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION PERCENT POINT FUNCTION .
C             VALUE PPF FOR THE CHI-SQUARED DISTRIBUTION
C             WITH DEGREES OF FREEDOM PARAMETER = NU.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
C     RESTRICTIONS--NU SHOULD BE A POSITIVE INTEGER VARIABLE.
C                 --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DEXP, DLOG.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS)
C               COMPARED TO THE KNOWN NU = 2 (EXPONENTIAL)
C               RESULTS, AGREEMENT WAS HAD OUT TO 6 SIGNIFICANT
C               DIGITS FOR ALL TESTED P IN THE RANGE P = .001 TO
C               P = .999.  FOR P = .95 AND SMALLER, THE AGREEMENT
C               WAS EVEN BETTER--7 SIGNIFICANT DIGITS.
C               (NOTE THAT THE TABULATED VALUES GIVEN IN THE WILK,
C               GNANADESIKAN, AND HUYETT REFERENCE BELOW, PAGE 20,
C               ARE IN ERROR FOR AT LEAST THE GAMMA = 1 CASE--
C               THE WORST DETECTED ERROR WAS AGREEMENT TO ONLY 3
C               SIGNIFICANT DIGITS (IN THEIR 8 SIGNIFICANT DIGIT TABLE)
C               FOR P = .999.)
C     REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY
C                 PLOTS FOR THE GAMMA DISTRIBUTION',
C                 TECHNOMETRICS, 1962, PAGES 1-15,
C                 ESPECIALLY PAGES 3-5.
C               --NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 257, FORMULA 6.1.41,
C                 AND PAGES 940-943.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--1, 1970, PAGES 166-206.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 46-51.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY.
C     LANGUAGE--ANSI FORTRAN (1966)
C     VERSION NUMBER--82/7
C     ORIGINAL VERSION--SEPTEMBER 1975.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --DECEMBER  1981.
C     UPDATED         --MAY       1982.
C     UPDATED         --JUNE      1987.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DP,DGAMMA
      DOUBLE PRECISION Z,Z2,DEN,A,B,C,D
      DOUBLE PRECISION XMIN0,XMIN,AI,XMAX,DX,PCALC,XMID
      DOUBLE PRECISION XLOWER,XUPPER,XDEL
      DOUBLE PRECISION SUM,TERM,CUT1,CUT2,AJ,CUTOFF,T
      DOUBLE PRECISION DLG,DLT,DLX,DLPCAL
      DOUBLE PRECISION DLP,DLGAMM,DLXMI0
      DOUBLE PRECISION Z2INV
      DOUBLE PRECISION DEXP,DLOG
C
      DIMENSION D(10)
C
C-----DATA STATEMENTS-------------------------------------------------
C
      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
C-----START POINT-----------------------------------------------------
C
      XMID=0.0
      XLOWER=0.0
      XUPPER=0.0
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      PPF=0.0
      IF(P.LT.0.0.OR.P.GE.1.0)THEN
        WRITE(IPR,1)
    1   FORMAT('***** ERROR IN CHSPPF--')
        WRITE(IPR,2)
    2   FORMAT('      THE FIRST ARGUMENT TO CHSPPF IS OUTSIDE ',
     1         'THE ALLOWABLE (0,1) INTERVAL')
        WRITE(IPR,46)P
   46   FORMAT('      THE VALUE OF THE ARGUMENT IS ',G15.7)
        GOTO9000
      ELSEIF(NU.LT.1)THEN
        WRITE(IPR,1)
        WRITE(IPR,15)
   15   FORMAT('      THE DEGREES OF FREEDOM ARGUMENT IS NON-POSITIVE')
        WRITE(IPR,47)NU
   47   FORMAT('      THE VALUE OF THE ARGUMENT IS ',I8)
        GOTO9000
      ENDIF
C
C     EXPRESS THE CHI-SQUARED DISTRIBUTION PERCENT POINT
C     FUNCTION IN TERMS OF THE EQUIVALENT GAMMA
C     DISTRIBUTION PERCENT POINT FUNCTION,
C     AND THEN EVALUATE THE LATTER.
C
      DP=P
      DNU=NU
      DGAMMA=DNU/2.0D0
      MAXIT=10000
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
      A=(Z-0.5D0)*DLOG(Z)-Z+C
      Z2INV=1.0D0/Z2
      B=D(9)
      B=Z2INV*B+D(8)
      B=Z2INV*B+D(7)
      B=Z2INV*B+D(6)
      B=Z2INV*B+D(5)
      B=Z2INV*B+D(4)
      B=Z2INV*B+D(3)
      B=Z2INV*B+D(2)
      B=Z2INV*B+D(1)
      B=(1.0D0/Z)*B
      DLG=(A+B)-DLOG(DEN)
C
C     DETERMINE LOWER AND UPPER LIMITS ON THE DESIRED 100P
C     PERCENT POINT.
C
      ILOOP=1
      DLP=DLOG(DP)
      DLGAMM=DLOG(DGAMMA)
      DLXMI0=(1.0D0/DGAMMA)*(DLP+DLGAMM+DLG)
      XMIN0=DEXP(DLXMI0)
      XMIN=XMIN0
      ICOUNT=1
  350 AI=ICOUNT
      XMAX=AI*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.0D0
C
C     NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED.
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.0D0
      GOTO590
  580 XUPPER=XMID
      XMID=(XMID+XLOWER)/2.0D0
  590 XDEL=XMID-XLOWER
      IF(XDEL.LT.0.0D0)XDEL=-XDEL
      ICOUNT=ICOUNT+1
      IF(XDEL.LT.0.0000000001D0.OR.ICOUNT.GT.100)GOTO570
      GOTO550
  570 PPF=2.0D0*XMID
      GOTO9000
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.0D0/DGAMMA
      TERM=1.0D0/DGAMMA
      CUT1=DX-DGAMMA
      CUT2=DX*10000000000.0D0
      DO700J=1,MAXIT
      AJ=J
      TERM=DX*TERM/(DGAMMA+AJ)
      SUM=SUM+TERM
      CUTOFF=CUT1+(CUT2*TERM/SUM)
      IF(AJ.GT.CUTOFF)GOTO750
  700 CONTINUE
      WRITE(IPR,1)
      WRITE(IPR,705)MAXIT
  705 FORMAT('     THE NUMBER OF ITERATIONS EXCEEDS ',I7)
      WRITE(IPR,706)P
  706 FORMAT('     THE INPUT VALUE OF P     IS ',G15.7)
      WRITE(IPR,707)NU
  707 FORMAT('     THE INPUT VALUE OF NU    IS ',I8)
      WRITE(IPR,708)
  708 FORMAT('     THE OUTPUT VALUE OF PPF HAS BEEN SET TO 0')
      PPF=0.0
      GOTO9000
C
  750 T=SUM
      DLT=DLOG(T)
      DLX=DLOG(DX)
      DLPCAL=DGAMMA*DLX-DX+DLT-DLG
      PCALC=DEXP(DLPCAL)
      IF(ILOOP.EQ.1)GOTO360
      GOTO560
C
 9000 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION D9LGMC (X, IPR)
C***BEGIN PROLOGUE  D9LGMC
C***SUBSIDIARY
C***PURPOSE  Compute the log Gamma correction factor so that
C            LOG(DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-5.)*LOG(X) - X
C            + D9LGMC(X).
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      DOUBLE PRECISION (R9LGMC-S, D9LGMC-D, C9LGMC-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, CORRECTION TERM, FNLIB,
C             LOG GAMMA, LOGARITHM, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute the log gamma correction factor for X .GE. 10. so that
C LOG (DGAMMA(X)) = LOG(SQRT(2*PI)) + (X-.5)*LOG(X) - X + D9lGMC(X)
C
C Series for ALGM       on the interval  0.          to  1.00000E-02
C                                        with weighted error   1.28E-31
C                                         log weighted error  30.89
C                               significant figures required  29.81
C                                    decimal places required  31.48
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900720  Routine changed from user-callable to subsidiary.  (WRB)
C***END PROLOGUE  D9LGMC
      DOUBLE PRECISION X, ALGMCS(15), XBIG, XMAX, DCSEVL
      LOGICAL FIRST
      SAVE ALGMCS, NALGM, XBIG, XMAX, FIRST
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA ALGMCS(  1) / +.1666389480 4518632472 0572965082 2 D+0      /
      DATA ALGMCS(  2) / -.1384948176 0675638407 3298605913 5 D-4      /
      DATA ALGMCS(  3) / +.9810825646 9247294261 5717154748 7 D-8      /
      DATA ALGMCS(  4) / -.1809129475 5724941942 6330626671 9 D-10     /
      DATA ALGMCS(  5) / +.6221098041 8926052271 2601554341 6 D-13     /
      DATA ALGMCS(  6) / -.3399615005 4177219443 0333059966 6 D-15     /
      DATA ALGMCS(  7) / +.2683181998 4826987489 5753884666 6 D-17     /
      DATA ALGMCS(  8) / -.2868042435 3346432841 4462239999 9 D-19     /
      DATA ALGMCS(  9) / +.3962837061 0464348036 7930666666 6 D-21     /
      DATA ALGMCS( 10) / -.6831888753 9857668701 1199999999 9 D-23     /
      DATA ALGMCS( 11) / +.1429227355 9424981475 7333333333 3 D-24     /
      DATA ALGMCS( 12) / -.3547598158 1010705471 9999999999 9 D-26     /
      DATA ALGMCS( 13) / +.1025680058 0104709120 0000000000 0 D-27     /
      DATA ALGMCS( 14) / -.3401102254 3167487999 9999999999 9 D-29     /
      DATA ALGMCS( 15) / +.1276642195 6300629333 3333333333 3 D-30     /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  D9LGMC
      IF (FIRST) THEN
         NALGM = INITDS (ALGMCS, 15, REAL(D1MACH(3)) )
         XBIG = 1.0D0/SQRT(D1MACH(3))
         XMAX = EXP (MIN(LOG(D1MACH(2)/12.D0), -LOG(12.D0*D1MACH(1))))
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LT. 10.D0) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM D9LGMC.  X MUST BE GREATER THAN ')
        WRITE(IPR,12)
   12   FORMAT('      OR EQUAL TO 10.                          ******')
        D9LGMC = 0.D0
        RETURN
      ENDIF
      IF (X.GE.XMAX) GO TO 20
C
      D9LGMC = 1.D0/(12.D0*X)
      IF (X.LT.XBIG) D9LGMC = DCSEVL (2.0D0*(10.D0/X)**2-1.D0, ALGMCS,
     1  NALGM, IPR) / X
      RETURN
C
 20   D9LGMC = 0.D0
      WRITE(IPR,21)
 21   FORMAT('***** WARNING FROM D9LGMC.  X SO BIG D9LCMC UNDERFLOWS.')
      RETURN
C
      END
      DOUBLE PRECISION FUNCTION DBETAI (X, PIN, QIN, IPR)
C***BEGIN PROLOGUE  DBETAI
C***PURPOSE  Calculate the incomplete Beta function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7F
C***TYPE      DOUBLE PRECISION (BETAI-S, DBETAI-D)
C***KEYWORDS  FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C   DBETAI calculates the DOUBLE PRECISION incomplete beta function.
C
C   The incomplete beta function ratio is the probability that a
C   random variable from a beta distribution having parameters PIN and
C   QIN will be less than or equal to X.
C
C     -- Input Arguments -- All arguments are DOUBLE PRECISION.
C   X      upper limit of integration.  X must be in (0,1) inclusive.
C   PIN    first beta distribution parameter.  PIN must be .GT. 0.0.
C   QIN    second beta distribution parameter.  QIN must be .GT. 0.0.
C
C***REFERENCES  Nancy E. Bosten and E. L. Battiste, Remark on Algorithm
C                 179, Communications of the ACM 17, 3 (March 1974),
C                 pp. 156.
C***ROUTINES CALLED  D1MACH, DLBETA, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
C***END PROLOGUE  DBETAI
      DOUBLE PRECISION X, PIN, QIN, ALNEPS, ALNSML, C, EPS, FINSUM, P,
     1  PS, Q, SML, TERM, XB, XI, Y, DLBETA, P1
      LOGICAL FIRST
      SAVE EPS, ALNEPS, SML, ALNSML, FIRST
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DBETAI
      IF (FIRST) THEN
         EPS = D1MACH(3)
         ALNEPS = LOG (EPS)
         SML = D1MACH(1)
         ALNSML = LOG (SML)
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LT. 0.D0 .OR. X .GT. 1.D0) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM DBETAI.  X IS NOT IN THE RANGE ')
        WRITE(IPR,12)
   12   FORMAT('      (0,1).                                    *****')
        DBETAI = 0.D0
        RETURN
      ENDIF
      IF (PIN .LE. 0.D0 .OR. QIN .LE. 0.D0) THEN
        WRITE(IPR,16)
   16   FORMAT('***** ERROR FROM DBETAI.  P AND/OR Q IS LESS THAN ')
        WRITE(IPR,17)
   17   FORMAT('      OR EQUAL TO ZERO.                         *****')
        DBETAI = 0.D0
        RETURN
      ENDIF
C
      Y = X
      P = PIN
      Q = QIN
      IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20
      IF (X.LT.0.2D0) GO TO 20
      Y = 1.0D0 - Y
      P = QIN
      Q = PIN
C
 20   IF ((P+Q)*Y/(P+1.D0).LT.EPS) GO TO 80
C
C EVALUATE THE INFINITE SUM FIRST.  TERM WILL EQUAL
C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) .
C
      PS = Q - AINT(Q)
      IF (PS.EQ.0.D0) PS = 1.0D0
      XB = P*LOG(Y) - DLBETA(PS,P,IPR) - LOG(P)
      DBETAI = 0.0D0
      IF (XB.LT.ALNSML) GO TO 40
C
      DBETAI = EXP (XB)
      TERM = DBETAI*P
      IF (PS.EQ.1.0D0) GO TO 40
      N = MAX (ALNEPS/LOG(Y), 4.0D0)
      DO 30 I=1,N
        XI = I
        TERM = TERM * (XI-PS)*Y/XI
        DBETAI = DBETAI + TERM/(P+XI)
 30   CONTINUE
C
C NOW EVALUATE THE FINITE SUM, MAYBE.
C
 40   IF (Q.LE.1.0D0) GO TO 70
C
      XB = P*LOG(Y) + Q*LOG(1.0D0-Y) - DLBETA(P,Q) - LOG(Q)
      IB = MAX (XB/ALNSML, 0.0D0)
      TERM = EXP(XB - IB*ALNSML)
      C = 1.0D0/(1.D0-Y)
      P1 = Q*C/(P+Q-1.D0)
C
      FINSUM = 0.0D0
      N = Q
      IF (Q.EQ.DBLE(N)) N = N - 1
      DO 50 I=1,N
        IF (P1.LE.1.0D0 .AND. TERM/EPS.LE.FINSUM) GO TO 60
        XI = I
        TERM = (Q-XI+1.0D0)*C*TERM/(P+Q-XI)
C
        IF (TERM.GT.1.0D0) IB = IB - 1
        IF (TERM.GT.1.0D0) TERM = TERM*SML
C
        IF (IB.EQ.0) FINSUM = FINSUM + TERM
 50   CONTINUE
C
 60   DBETAI = DBETAI + FINSUM
 70   IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI
      DBETAI = MAX (MIN (DBETAI, 1.0D0), 0.0D0)
      RETURN
C
 80   DBETAI = 0.0D0
      XB = P*LOG(MAX(Y,SML)) - LOG(P) - DLBETA(P,Q,IPR)
      IF (XB.GT.ALNSML .AND. Y.NE.0.0D0) DBETAI = EXP(XB)
      IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION DCSEVL (X, CS, N, IPR)
C***BEGIN PROLOGUE  DCSEVL
C***PURPOSE  Evaluate a Chebyshev series.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C3A2
C***TYPE      DOUBLE PRECISION (CSEVL-S, DCSEVL-D)
C***KEYWORDS  CHEBYSHEV SERIES, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C  Evaluate the N-term Chebyshev series CS at X.  Adapted from
C  a method presented in the paper by Broucke referenced below.
C
C       Input Arguments --
C  X    value at which the series is to be evaluated.
C  CS   array of N terms of a Chebyshev series.  In evaluating
C       CS, only half the first coefficient is summed.
C  N    number of terms in array CS.
C
C***REFERENCES  R. Broucke, Ten subroutines for the manipulation of
C                 Chebyshev series, Algorithm 446, Communications of
C                 the A.C.M. 16, (1973) pp. 254-256.
C               L. Fox and I. B. Parker, Chebyshev Polynomials in
C                 Numerical Analysis, Oxford University Press, 1968,
C                 page 56.
C***ROUTINES CALLED  D1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900329  Prologued revised extensively and code rewritten to allow
C           X to be slightly outside interval (-1,+1).  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DCSEVL
      DOUBLE PRECISION B0, B1, B2, CS(*), ONEPL, TWOX, X
      LOGICAL FIRST
      SAVE FIRST, ONEPL
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DCSEVL
      IF (FIRST) ONEPL = 1.0D0 + D1MACH(4)
      FIRST = .FALSE.
      IF (N .LT. 1) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM DCSEVL.  THE NUMBER OF TERMS IS ')
        WRITE(IPR,12)
   12   FORMAT('      LESS THAN OR EQUAL TO ZERO.                *****')
        DCSEVL = 0.D0
        RETURN
      ENDIF
      IF (N .GT. 1000) THEN
        WRITE(IPR,21)
   21   FORMAT('***** ERROR FROM DCSEVL.  THE NUMBER OF TERMS IS ')
        WRITE(IPR,22)
   22   FORMAT('      GREATER THAN 1000.                         *****')
        DCSEVL = 0.D0
        RETURN
      ENDIF
      IF (ABS(X) .GT. ONEPL) THEN
        WRITE(IPR,31)
   31   FORMAT('***** WARNING FROM DCSEVL.  X IS OUTSIDE THE ')
        WRITE(IPR,32)
   32   FORMAT('      INTERVAL (-1,+1).                          *****')
      ENDIF
C
      B1 = 0.0D0
      B0 = 0.0D0
      TWOX = 2.0D0*X
      DO 10 I = 1,N
         B2 = B1
         B1 = B0
         NI = N + 1 - I
         B0 = TWOX*B1 - B2 + CS(NI)
   10 CONTINUE
C
      DCSEVL = 0.5D0*(B0-B2)
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION DLBETA (A, B, IPR)
C***BEGIN PROLOGUE  DLBETA
C***PURPOSE  Compute the natural logarithm of the complete Beta
C            function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7B
C***TYPE      DOUBLE PRECISION (ALBETA-S, DLBETA-D, CLBETA-C)
C***KEYWORDS  FNLIB, LOGARITHM OF THE COMPLETE BETA FUNCTION,
C             SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DLBETA(A,B) calculates the double precision natural logarithm of
C the complete beta function for double precision arguments
C A and B.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D9LGMC, DGAMMA, DLNGAM, DLNREL, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900727  Added EXTERNAL statement.  (WRB)
C***END PROLOGUE  DLBETA
      DOUBLE PRECISION A, B, P, Q, CORR, SQ2PIL, D9LGMC, DGAMMA, DLNGAM,
     1  DLNREL
      EXTERNAL DGAMMA
      SAVE SQ2PIL
C
      DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
C***FIRST EXECUTABLE STATEMENT  DLBETA
      P = MIN (A, B)
      Q = MAX (A, B)
C
      IF (P .LE. 0.D0) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM DLBETA.  BOTH INPUT ARGUMENTS ')
        WRITE(IPR,12)
   12   FORMAT('      MUST BE GREATER THAN ZERO.               ******')
        DLBETA = 0.D0
        RETURN
      ENDIF
C
      IF (P.GE.10.D0) GO TO 30
      IF (Q.GE.10.D0) GO TO 20
C
C P AND Q ARE SMALL.
C
      DLBETA = LOG (DGAMMA(P,IPR) * (DGAMMA(Q,IPR)/DGAMMA(P+Q,IPR)) )
      RETURN
C
C P IS SMALL, BUT Q IS BIG.
C
 20   CORR = D9LGMC(Q,IPR) - D9LGMC(P+Q,IPR)
      DLBETA = DLNGAM(P,IPR) + CORR + P - P*LOG(P+Q)
     1  + (Q-0.5D0)*DLNREL(-P/(P+Q),IPR)
      RETURN
C
C P AND Q ARE BIG.
C
 30   CORR = D9LGMC(P,IPR) + D9LGMC(Q,IPR) - D9LGMC(P+Q,IPR)
      DLBETA = -0.5D0*LOG(Q) + SQ2PIL + CORR + (P-0.5D0)*LOG(P/(P+Q))
     1  + Q*DLNREL(-P/(P+Q),IPR)
      RETURN
C
      END
      DOUBLE PRECISION FUNCTION DLNGAM (X, IPR)
C***BEGIN PROLOGUE  DLNGAM
C***PURPOSE  Compute the logarithm of the absolute value of the Gamma
C            function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A
C***TYPE      DOUBLE PRECISION (ALNGAM-S, DLNGAM-D, CLNGAM-C)
C***KEYWORDS  ABSOLUTE VALUE, COMPLETE GAMMA FUNCTION, FNLIB, LOGARITHM,
C             SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DLNGAM(X) calculates the double precision logarithm of the
C absolute value of the Gamma function for double precision
C argument X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, D9LGMC, DGAMMA, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900727  Added EXTERNAL statement.  (WRB)
C***END PROLOGUE  DLNGAM
      DOUBLE PRECISION X, DXREL, PI, SINPIY, SQPI2L, SQ2PIL, XMAX,
     1  Y, DGAMMA, D9LGMC, TEMP
      LOGICAL FIRST
      EXTERNAL DGAMMA
      SAVE SQ2PIL, SQPI2L, PI, XMAX, DXREL, FIRST
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
      DATA SQPI2L / +.2257913526 4472743236 3097614947 441 D+0    /
      DATA PI / 3.1415926535 8979323846 2643383279 50 D0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DLNGAM
      IF (FIRST) THEN
         TEMP = 1.D0/LOG(D1MACH(2))
         XMAX = TEMP*D1MACH(2)
         DXREL = SQRT(D1MACH(4))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS (X)
      IF (Y.GT.10.D0) GO TO 20
C
C LOG (ABS (DGAMMA(X)) ) FOR ABS(X) .LE. 10.0
C
      DLNGAM = LOG (ABS (DGAMMA(X, IPR)) )
      RETURN
C
C LOG ( ABS (DGAMMA(X)) ) FOR ABS(X) .GT. 10.0
C
 20   IF (Y .GT. XMAX) THEN
        WRITE(IPR,21)
   21   FORMAT('***** ERROR FROM DLNGAM.  ABSOLUTE VALUE OF X SO ')
        WRITE(IPR,22)
   22   FORMAT('      LARGE THAT DLNGAM OVERFLOWS.             ******')
        DLNGAM = 0.D0
        RETURN
      ENDIF
C
      IF (X.GT.0.D0) DLNGAM = SQ2PIL + (X-0.5D0)*LOG(X) - X +
     1    D9LGMC(Y,IPR)
      IF (X.GT.0.D0) RETURN
C
      SINPIY = ABS (SIN(PI*Y))
      IF (SINPIY .EQ. 0.D0) THEN
        WRITE(IPR,31)
   31   FORMAT('***** ERROR FROM DLNGAM.  X IS A NEGATIVE INTEGER. ')
        DLNGAM = 0.D0
        RETURN
      ENDIF
C
      IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL)THEN
        WRITE(IPR,41)
   41   FORMAT('***** WARNING FROM DLNGAM.  ANSWER LESS THAN HALF ')
        WRITE(IPR,42)
   42   FORMAT('      PRECISION BECAUSE X IS TOO NEAR A NEGATIVE ')
        WRITE(IPR,43)
   43   FORMAT('      INTEGER.                                 *****')
      ENDIF
C
      DLNGAM = SQPI2L + (X-0.5D0)*LOG(Y) - X - LOG(SINPIY) -
     1         D9LGMC(Y, IPR)
      RETURN
C
      END
      DOUBLE PRECISION FUNCTION DLNREL (X, IPR)
C***BEGIN PROLOGUE  DLNREL
C***PURPOSE  Evaluate ln(1+X) accurate in the sense of relative error.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4B
C***TYPE      DOUBLE PRECISION (ALNREL-S, DLNREL-D, CLNREL-C)
C***KEYWORDS  ELEMENTARY FUNCTIONS, FNLIB, LOGARITHM
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DLNREL(X) calculates the double precision natural logarithm of
C (1.0+X) for double precision argument X.  This routine should
C be used when X is small and accurate to calculate the logarithm
C accurately (in the relative error sense) in the neighborhood
C of 1.0.
C
C Series for ALNR       on the interval -3.75000E-01 to  3.75000E-01
C                                        with weighted error   6.35E-32
C                                         log weighted error  31.20
C                               significant figures required  30.93
C                                    decimal places required  32.01
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DLNREL
      DOUBLE PRECISION ALNRCS(43), X, XMIN,  DCSEVL
      LOGICAL FIRST
      SAVE ALNRCS, NLNREL, XMIN, FIRST
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA ALNRCS(  1) / +.1037869356 2743769800 6862677190 98 D+1     /
      DATA ALNRCS(  2) / -.1336430150 4908918098 7660415531 33 D+0     /
      DATA ALNRCS(  3) / +.1940824913 5520563357 9261993747 50 D-1     /
      DATA ALNRCS(  4) / -.3010755112 7535777690 3765377765 92 D-2     /
      DATA ALNRCS(  5) / +.4869461479 7154850090 4563665091 37 D-3     /
      DATA ALNRCS(  6) / -.8105488189 3175356066 8099430086 22 D-4     /
      DATA ALNRCS(  7) / +.1377884779 9559524782 9382514960 59 D-4     /
      DATA ALNRCS(  8) / -.2380221089 4358970251 3699929149 35 D-5     /
      DATA ALNRCS(  9) / +.4164041621 3865183476 3918599019 89 D-6     /
      DATA ALNRCS( 10) / -.7359582837 8075994984 2668370319 98 D-7     /
      DATA ALNRCS( 11) / +.1311761187 6241674949 1522943450 11 D-7     /
      DATA ALNRCS( 12) / -.2354670931 7742425136 6960923301 75 D-8     /
      DATA ALNRCS( 13) / +.4252277327 6034997775 6380529625 67 D-9     /
      DATA ALNRCS( 14) / -.7719089413 4840796826 1081074933 00 D-10    /
      DATA ALNRCS( 15) / +.1407574648 1359069909 2153564721 91 D-10    /
      DATA ALNRCS( 16) / -.2576907205 8024680627 5370786275 84 D-11    /
      DATA ALNRCS( 17) / +.4734240666 6294421849 1543950059 38 D-12    /
      DATA ALNRCS( 18) / -.8724901267 4742641745 3012632926 75 D-13    /
      DATA ALNRCS( 19) / +.1612461490 2740551465 7398331191 15 D-13    /
      DATA ALNRCS( 20) / -.2987565201 5665773006 7107924168 15 D-14    /
      DATA ALNRCS( 21) / +.5548070120 9082887983 0413216972 79 D-15    /
      DATA ALNRCS( 22) / -.1032461915 8271569595 1413339619 32 D-15    /
      DATA ALNRCS( 23) / +.1925023920 3049851177 8785032448 68 D-16    /
      DATA ALNRCS( 24) / -.3595507346 5265150011 1897078442 66 D-17    /
      DATA ALNRCS( 25) / +.6726454253 7876857892 1945742267 73 D-18    /
      DATA ALNRCS( 26) / -.1260262416 8735219252 0824256375 46 D-18    /
      DATA ALNRCS( 27) / +.2364488440 8606210044 9161589555 19 D-19    /
      DATA ALNRCS( 28) / -.4441937705 0807936898 8783891797 33 D-20    /
      DATA ALNRCS( 29) / +.8354659446 4034259016 2412939946 66 D-21    /
      DATA ALNRCS( 30) / -.1573155941 6479562574 8992535210 66 D-21    /
      DATA ALNRCS( 31) / +.2965312874 0247422686 1543697066 66 D-22    /
      DATA ALNRCS( 32) / -.5594958348 1815947292 1560132266 66 D-23    /
      DATA ALNRCS( 33) / +.1056635426 8835681048 1872841386 66 D-23    /
      DATA ALNRCS( 34) / -.1997248368 0670204548 3149994666 66 D-24    /
      DATA ALNRCS( 35) / +.3778297781 8839361421 0498559999 99 D-25    /
      DATA ALNRCS( 36) / -.7153158688 9081740345 0381653333 33 D-26    /
      DATA ALNRCS( 37) / +.1355248846 3674213646 5020245333 33 D-26    /
      DATA ALNRCS( 38) / -.2569467304 8487567430 0798293333 33 D-27    /
      DATA ALNRCS( 39) / +.4874775606 6216949076 4595199999 99 D-28    /
      DATA ALNRCS( 40) / -.9254211253 0849715321 1323733333 33 D-29    /
      DATA ALNRCS( 41) / +.1757859784 1760239233 2697600000 00 D-29    /
      DATA ALNRCS( 42) / -.3341002667 7731010351 3770666666 66 D-30    /
      DATA ALNRCS( 43) / +.6353393618 0236187354 1802666666 66 D-31    /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DLNREL
      IF (FIRST) THEN
         NLNREL = INITDS (ALNRCS, 43, 0.1*REAL(D1MACH(3)), IPR)
         XMIN = -1.0D0 + SQRT(D1MACH(4))
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LE. (-1.D0)) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM DLNREL.  X IS LESS THAN OR ')
        WRITE(IPR,12)
   12   FORMAT('      EQUAL TO -1.                             ******')
        DLNREL = 0.0
        RETURN
      ENDIF
      IF (X .LT. XMIN) THEN
        WRITE(IPR,21)
 21     FORMAT('***** WARNING FROM DLNREL.  ANSWER LESS THAN HALF ')
        WRITE(IPR,22)
 22     FORMAT('      PRECISION BECAUSE X IS TOO NEAR -1.     *****')
      ENDIF
C
      IF (ABS(X).LE.0.375D0) DLNREL = X*(1.D0 -
     1  X*DCSEVL (X/.375D0, ALNRCS, NLNREL, IPR))
C
      IF (ABS(X).GT.0.375D0) DLNREL = LOG (1.0D0+X)
C
      RETURN
      END
      DOUBLE PRECISION FUNCTION DGAMMA (X,IPR)
C***BEGIN PROLOGUE  DGAMMA
C***PURPOSE  Compute the complete Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A
C***TYPE      DOUBLE PRECISION (GAMMA-S, DGAMMA-D, CGAMMA-C)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DGAMMA(X) calculates the double precision complete Gamma function
C for double precision argument X.
C
C Series for GAM        on the interval  0.          to  1.00000E+00
C                                        with weighted error   5.79E-32
C                                         log weighted error  31.24
C                               significant figures required  30.00
C                                    decimal places required  32.05
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, D9LGMC, DCSEVL, DGAMLM, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920618  Removed space from variable name.  (RWC, WRB)
C***END PROLOGUE  DGAMMA
      DOUBLE PRECISION X, GAMCS(42), DXREL, PI, SINPIY, SQ2PIL, XMAX,
     1  XMIN, Y, D9LGMC, DCSEVL
      LOGICAL FIRST
C
      SAVE GAMCS, PI, SQ2PIL, NGAM, XMIN, XMAX, DXREL, FIRST
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
      DATA GAMCS(  1) / +.8571195590 9893314219 2006239994 2 D-2      /
      DATA GAMCS(  2) / +.4415381324 8410067571 9131577165 2 D-2      /
      DATA GAMCS(  3) / +.5685043681 5993633786 3266458878 9 D-1      /
      DATA GAMCS(  4) / -.4219835396 4185605010 1250018662 4 D-2      /
      DATA GAMCS(  5) / +.1326808181 2124602205 8400679635 2 D-2      /
      DATA GAMCS(  6) / -.1893024529 7988804325 2394702388 6 D-3      /
      DATA GAMCS(  7) / +.3606925327 4412452565 7808221722 5 D-4      /
      DATA GAMCS(  8) / -.6056761904 4608642184 8554829036 5 D-5      /
      DATA GAMCS(  9) / +.1055829546 3022833447 3182350909 3 D-5      /
      DATA GAMCS( 10) / -.1811967365 5423840482 9185589116 6 D-6      /
      DATA GAMCS( 11) / +.3117724964 7153222777 9025459316 9 D-7      /
      DATA GAMCS( 12) / -.5354219639 0196871408 7408102434 7 D-8      /
      DATA GAMCS( 13) / +.9193275519 8595889468 8778682594 0 D-9      /
      DATA GAMCS( 14) / -.1577941280 2883397617 6742327395 3 D-9      /
      DATA GAMCS( 15) / +.2707980622 9349545432 6654043308 9 D-10     /
      DATA GAMCS( 16) / -.4646818653 8257301440 8166105893 3 D-11     /
      DATA GAMCS( 17) / +.7973350192 0074196564 6076717535 9 D-12     /
      DATA GAMCS( 18) / -.1368078209 8309160257 9949917230 9 D-12     /
      DATA GAMCS( 19) / +.2347319486 5638006572 3347177168 8 D-13     /
      DATA GAMCS( 20) / -.4027432614 9490669327 6657053469 9 D-14     /
      DATA GAMCS( 21) / +.6910051747 3721009121 3833697525 7 D-15     /
      DATA GAMCS( 22) / -.1185584500 2219929070 5238712619 2 D-15     /
      DATA GAMCS( 23) / +.2034148542 4963739552 0102605193 2 D-16     /
      DATA GAMCS( 24) / -.3490054341 7174058492 7401294910 8 D-17     /
      DATA GAMCS( 25) / +.5987993856 4853055671 3505106602 6 D-18     /
      DATA GAMCS( 26) / -.1027378057 8722280744 9006977843 1 D-18     /
      DATA GAMCS( 27) / +.1762702816 0605298249 4275966074 8 D-19     /
      DATA GAMCS( 28) / -.3024320653 7353062609 5877211204 2 D-20     /
      DATA GAMCS( 29) / +.5188914660 2183978397 1783355050 6 D-21     /
      DATA GAMCS( 30) / -.8902770842 4565766924 4925160106 6 D-22     /
      DATA GAMCS( 31) / +.1527474068 4933426022 7459689130 6 D-22     /
      DATA GAMCS( 32) / -.2620731256 1873629002 5732833279 9 D-23     /
      DATA GAMCS( 33) / +.4496464047 8305386703 3104657066 6 D-24     /
      DATA GAMCS( 34) / -.7714712731 3368779117 0390152533 3 D-25     /
      DATA GAMCS( 35) / +.1323635453 1260440364 8657271466 6 D-25     /
      DATA GAMCS( 36) / -.2270999412 9429288167 0231381333 3 D-26     /
      DATA GAMCS( 37) / +.3896418998 0039914493 2081663999 9 D-27     /
      DATA GAMCS( 38) / -.6685198115 1259533277 9212799999 9 D-28     /
      DATA GAMCS( 39) / +.1146998663 1400243843 4761386666 6 D-28     /
      DATA GAMCS( 40) / -.1967938586 3451346772 9510399999 9 D-29     /
      DATA GAMCS( 41) / +.3376448816 5853380903 3489066666 6 D-30     /
      DATA GAMCS( 42) / -.5793070335 7821357846 2549333333 3 D-31     /
      DATA PI / 3.1415926535 8979323846 2643383279 50 D0 /
      DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DGAMMA
      IF (FIRST) THEN
         NGAM = INITDS (GAMCS, 42, 0.1*REAL(D1MACH(3)), IPR )
C
         CALL DGAMLM (XMIN, XMAX, IPR)
         DXREL = SQRT(D1MACH(4))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS(X)
      IF (Y.GT.10.D0) GO TO 50
C
C COMPUTE GAMMA(X) FOR -XBND .LE. X .LE. XBND.  REDUCE INTERVAL AND FIND
C GAMMA(1+Y) FOR 0.0 .LE. Y .LT. 1.0 FIRST OF ALL.
C
      N = X
      IF (X.LT.0.D0) N = N - 1
      Y = X - N
      N = N - 1
      DGAMMA = 0.9375D0 + DCSEVL (2.D0*Y-1.D0, GAMCS, NGAM, IPR)
      IF (N.EQ.0) RETURN
C
      IF (N.GT.0) GO TO 30
C
C COMPUTE GAMMA(X) FOR X .LT. 1.0
C
      N = -N
      IF (X .EQ. 0.D0) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM DGAMMA.  X IS 0.  ******')
        DGAMMA = 0.D0
        RETURN
      ENDIF
      IF (X .LT. 0.0 .AND. X+N-2 .EQ. 0.D0)THEN
        WRITE(IPR,16)
   16   FORMAT('***** ERROR FROM DGAMMA.  X IS A NEGATIVE INTEGER. ',
     1         '****')
        DGAMMA = 0.D0
        RETURN
      ENDIF
      IF(X .LT. (-0.5D0) .AND. ABS((X-AINT(X-0.5D0))/X) .LT. DXREL)THEN
        WRITE(IPR,21)
   21   FORMAT('***** WARNING FROM DGAMMA.  ANSWER IS LESS THAN ')
      ENDIF
   22 FORMAT('      HALF PRECISION BECAUSE X IS TOO NEAR A ')
   23 FORMAT('      NEGATIVE INTEGER.                          *****')
C
      DO 20 I=1,N
        DGAMMA = DGAMMA/(X+I-1 )
 20   CONTINUE
      RETURN
C
C GAMMA(X) FOR X .GE. 2.0 AND X .LE. 10.0
C
 30   DO 40 I=1,N
        DGAMMA = (Y+I) * DGAMMA
 40   CONTINUE
      RETURN
C
C GAMMA(X) FOR ABS(X) .GT. 10.0.  RECALL Y = ABS(X).
C
 50   IF (X .GT. XMAX) THEN
        WRITE(IPR,51)
   51   FORMAT('***** ERROR FROM DGAMMA.  X IS SO BIG THAT THE ')
        WRITE(IPR,52)
   52   FORMAT('      DGAMMA FUNCTION OVERFLOWS.               *****')
        DGAMMA = 0.D0
        RETURN
      ENDIF
C
      DGAMMA = 0.D0
      IF (X .LT. XMIN) THEN
        WRITE(IPR,56)
   56   FORMAT('***** WARNING FROM DGAMMA.  X IS SO SMALL THAT THE ')
        WRITE(IPR,57)
   57   FORMAT('      DGAMMA FUNCTION UNDERFLOWS.              *****')
      ENDIF
      IF (X.LT.XMIN) RETURN
C
      DGAMMA = EXP ((Y-0.5D0)*LOG(Y) - Y + SQ2PIL + D9LGMC(Y, IPR) )
      IF (X.GT.0.D0) RETURN
C
      IF (ABS((X-AINT(X-0.5D0))/X) .LT. DXREL) THEN
        WRITE(IPR,61)
   61   FORMAT('***** WARNING FROM DGAMMA.  ANSWER IS LESS THAN ')
        WRITE(IPR,62)
   62   FORMAT('      PRECISION BECAUSE X IS TOO NEAR A NEGATIVE ')
        WRITE(IPR,63)
   63   FORMAT('      NUMBER.                                 *****')
      ENDIF
C
      SINPIY = SIN (PI*Y)
      IF (SINPIY .EQ. 0.D0) THEN
        WRITE(IPR,71)
   71   FORMAT('***** ERROR FROM DGAMMA.  X IS A NEGATIVE INTEGER. ',
     1         '****')
        DGAMMA = 0.D0
        RETURN
      ENDIF
C
      DGAMMA = -PI/(Y*SINPIY*DGAMMA)
C
      RETURN
      END
      SUBROUTINE DGAMLM (XMIN, XMAX, IPR)
C***BEGIN PROLOGUE  DGAMLM
C***PURPOSE  Compute the minimum and maximum bounds for the argument in
C            the Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A, R2
C***TYPE      DOUBLE PRECISION (GAMLIM-S, DGAMLM-D)
C***KEYWORDS  COMPLETE GAMMA FUNCTION, FNLIB, LIMITS, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Calculate the minimum and maximum legal bounds for X in gamma(X).
C XMIN and XMAX are not the only bounds, but they are the only non-
C trivial ones to calculate.
C
C             Output Arguments --
C XMIN   double precision minimum legal value of X in gamma(X).  Any
C        smaller value of X might result in underflow.
C XMAX   double precision maximum legal value of X in gamma(X).  Any
C        larger value of X might cause overflow.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DGAMLM
      DOUBLE PRECISION XMIN, XMAX, ALNBIG, ALNSML, XLN, XOLD
C
C-----COMMON----------------------------------------------------------
C
C-----THIS IS COMMON/INCLUDE FILE DPCOMC.INC-----
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
C***FIRST EXECUTABLE STATEMENT  DGAMLM
      ALNSML = LOG(D1MACH(1))
      XMIN = -ALNSML
      DO 10 I=1,10
        XOLD = XMIN
        XLN = LOG(XMIN)
        XMIN = XMIN - XMIN*((XMIN+0.5D0)*XLN - XMIN - 0.2258D0 + ALNSML)
     1    / (XMIN*XLN+0.5D0)
        IF (ABS(XMIN-XOLD).LT.0.005D0) GO TO 20
 10   CONTINUE
      WRITE(IPR,11)
 11   FORMAT('***** ERROR FROM DGAMLM.  UNABLE TO FIND XMIN.  ******')
      RETURN
C
 20   XMIN = -XMIN + 0.01D0
C
      ALNBIG = LOG (D1MACH(2))
      XMAX = ALNBIG
      DO 30 I=1,10
        XOLD = XMAX
        XLN = LOG(XMAX)
        XMAX = XMAX - XMAX*((XMAX-0.5D0)*XLN - XMAX + 0.9189D0 - ALNBIG)
     1    / (XMAX*XLN-0.5D0)
        IF (ABS(XMAX-XOLD).LT.0.005D0) GO TO 40
 30   CONTINUE
      WRITE(IPR,21)
 21   FORMAT('***** ERROR FROM DGAMLM.  UNABLE TO FIND XMAX.  ******')
      RETURN
C
 40   XMAX = XMAX - 0.01D0
      XMIN = MAX (XMIN, -XMAX+1.D0)
C
      RETURN
      END
      SUBROUTINE DISTIN(X,NX,Y,NY,IPR,IDEBUG,IERROR)
C
C     PURPOSE--COMPUTE DISTINCT VALUES OF A VARIABLE--
C              Y(1) = X(1)
C              Y(2) = X(2) OR X(3) OR X(4) ETC., THE FIRST ONE
C                     OF WHICH IS DIFFERENT FROM Y(1);
C              Y(3) = X(3) OR X(4) OR X(5) ETC., THE FIRST ONE
C                     OF WHICH IS DIFFERENT FROM Y(1) AND Y(2);
C              ETC.
C     NOTE--IT IS PERMISSIBLE TO HAVE THE OUTPUT VECTOR Y(.)
C           BEING IDENTICAL TO THE INPUT VECTOR X(.).
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82/7
C     ORIGINAL VERSION--FEBRUARY  1979.
C     UPDATED         --APRIL     1979.
C     UPDATED         --AUGUST    1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
      DIMENSION Y(*)
C
C---------------------------------------------------------------------
C
C-----START POINT-----------------------------------------------------
C
      IERROR='NO'
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF DISTIN--')
        WRITE(IPR,52)IDEBUG,NX
   52   FORMAT('IDEBUG,NX = ',A4,2X,I8)
        DO55I=1,NX
          WRITE(IPR,56)I,X(I)
   56     FORMAT('I,X(I) = ',I8,E15.7)
   55   CONTINUE
      ENDIF
C
C               ********************************
C               **  COMPUTE DISTINCT VALUES.  **
C               ********************************
C
      NY=0
      IF(NX.LT.1)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,151)
  151   FORMAT('***** ERROR IN DISTIN--')
        WRITE(IPR,152)
  152   FORMAT('      THE NUMBER OF OBSERVATIONS FOR THE RESPONSE')
        WRITE(IPR,153)
  153   FORMAT('      VARIABLE IS LESS THAN ONE.')
        WRITE(IPR,157)NX
  157   FORMAT('      THE NUMBER OF OBSERVATIONS HERE = ',I8,'.')
      ELSE
        NY=1
        Y(NY)=X(1)
        IF(NX.LT.2)GOTO9000
        DO100I=2,NX
          DO120J=1,NY
            IF(X(I).EQ.Y(J))GOTO100
  120     CONTINUE
          NY=NY+1
          Y(NY)=X(I)
  100   CONTINUE
      ENDIF
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT.      **
C               *****************
C
 9000 CONTINUE
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF DISTIN--')
        WRITE(IPR,9012)IDEBUG,IERROR
 9012   FORMAT('IDEBUG,IERROR = ',A4,2X,A4)
        WRITE(IPR,9013)NX,NY
 9013   FORMAT('NX,NY = ',2I8)
        DO9015I=1,NY
          WRITE(IPR,9016)I,Y(I)
 9016     FORMAT('I,Y(I) = ',I8,G15.7)
 9015   CONTINUE
      ENDIF
C
      RETURN
      END
      FUNCTION INITDS (OS, NOS, ETA, IPR)
C***BEGIN PROLOGUE  INITDS
C***PURPOSE  Determine the number of terms needed in an orthogonal
C            polynomial series so that it meets a specified accuracy.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C3A2
C***TYPE      DOUBLE PRECISION (INITS-S, INITDS-D)
C***KEYWORDS  CHEBYSHEV, FNLIB, INITIALIZE, ORTHOGONAL POLYNOMIAL,
C             ORTHOGONAL SERIES, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C  Initialize the orthogonal series, represented by the array OS, so
C  that INITDS is the number of terms needed to insure the error is no
C  larger than ETA.  Ordinarily, ETA will be chosen to be one-tenth
C  machine precision.
C
C             Input Arguments --
C   OS     double precision array of NOS coefficients in an orthogonal
C          series.
C   NOS    number of coefficients in OS.
C   ETA    single precision scalar containing requested accuracy of
C          series.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891115  Modified error message.  (WRB)
C   891115  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  INITDS
      DOUBLE PRECISION OS(*)
C
C---------------------------------------------------------------------
C
C***FIRST EXECUTABLE STATEMENT  INITDS
      IF (NOS .LT. 1) THEN
        WRITE(IPR,11)
   11   FORMAT('***** ERROR FROM INITDS.  THE NUMBER OF ')
        WRITE(IPR,12)
   12   FORMAT('      COEFFICIENTS IS LESS THAN 1.      *****')
        INITDS = 0.D0
        RETURN
      ENDIF
C
      ERR = 0.
      DO 10 II = 1,NOS
        I = NOS + 1 - II
        ERR = ERR + ABS(REAL(OS(I)))
        IF (ERR.GT.ETA) GO TO 20
   10 CONTINUE
C
   20 IF (I .EQ. NOS) THEN
        WRITE(IPR,21)
 21     FORMAT('***** ERROR FROM INITDS.  CHEBYSHEV SERIES TOO ')
        WRITE(IPR,22)
 22     FORMAT('      SHORT FOR SPECIFIED ACCURACY.           *****')
      ENDIF
      INITDS = I
C
      RETURN
      END
      SUBROUTINE INITMC(IDEBUG,IPR)
C
C     PURPOSE--DEFINE MACHINE CONSTANTS (INTEGER, REAL, AND
C              DOUBLE PRECISION) FOR A PARTICULAR COMPUTER,
C
C              FOR PURPOSES OF LIMITS OF DETECTION, STRIP OUT
C              CONSTANTS FOR OLDER, OBSOLETE MACHINES
C
C     NOTE--THIS SUBROUTINE DOES NOT ADHERE TO 1966 ANSI STANDARD
C           OR THE 1977 ANSI STANDARD
C           BECAUSE IT USES OCTAL CONSTANTS.
C     NOTE--TO ALTER THIS SUBROUTINE FOR A PARTICULAR ENVIRONMENT,
C     THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C     REMOVING THE C FROM COLUMN 1.
C     NOTE--FOR IMPLEMENTATION CONVIENENCE, THE COMMENT LINES THAT
C     NEED TO BE DEACTIVATED ARE CODED AS "CXXXX" WHERE XXXX DEFINES
C     A PARTICULAR MACHINE.  THIS MEANS A SINGLE GLOBAL REPLACE CAN
C     BE USED TO UNCOMMENT THE APPROPRIATE LINES FOR A PARTICULAR
C     MACHINE (E.G., CHANGE '     ' TO '    ').
C
C               **************************************************
C               **  DESCRIPTION OF INTEGER MACHINE CONSTANTS    **
C               **************************************************
C
C     TO DESCRIBE I/O UNIT NUMBERS--
C
C       I1MACH( 1) = THE STANDARD INPUT UNIT.
C       I1MACH( 2) = THE STANDARD OUTPUT UNIT.
C       I1MACH( 3) = THE STANDARD PUNCH UNIT.
C       I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
C
C     TO DESCRIBE WORDS--
C
C       I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
C       I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT.
C
C     TO DESCRIBE INTEGERS--
C
C       ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
C                  SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C                  WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
C
C       I1MACH( 7) = A, THE BASE.
C       I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
C       I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
C
C     TO DESCIBE FLOATING-POINT NUMBERS--
C
C       ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C       BASE-B FORM
C                  SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C                  WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
C                  0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
C
C       I1MACH(10) = B, THE BASE.
C
C     TO DESCIBE SINGLE-PRECISION--
C
C       I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
C       I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
C       I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
C
C     TO DESCRIBE DOUBLE-PRECISION--
C
C       I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
C       I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
C       I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
C
C     THE VALUES OF
C     I1MACH(1) TO I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C     WITH THE LOCAL OPERATING SYSTEM.
C
C               *************************************************************
C               **  DESCRIPTION OF REAL (FLOATING POINT) MACHINE CONSTANTS  *
C               *************************************************************
C
C     R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C     R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C     R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C     R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C     R1MACH(5) = LOG10(B)
C
C     WHERE POSSIBLE, OCTAL OR HEXADECIMAL CONSTANTS HAVE BEEN USED
C     TO SPECIFY THE CONSTANTS EXACTLY WHICH HAS IN SOME CASES
C     REQUIRED THE USE OF EQUIVALENT INTEGER ARRAYS.
C
C               *********************************************************
C               **  DESCRIPTION OF DOUBLE PRECISION MACHINE CONSTANTS  **
C               *********************************************************
C
C     D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C     D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C     D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C     D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C     D1MACH( 5) = LOG10(B)
C
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE--301-921-3651
C     NOTE--THIS SUBROUTINE IS A MODIFICATION OF CODE
C           PROVIDED IN THE FOLLOWING ARTICLE--
C           CACM, 19XX.
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--SEPTEMBER 1980
C     UPDATED         --JULY      1981.
C     UPDATED         --AUGUST    1981.
C     UPDATED         --MAY       1982.
C     UPDATED         --JULY      1986 (IBM-PC AND MACINTOSH)
C     UPDATED         --OCTOBER   1986 (SUN)
C     UPDATED         --FEBRUARY  1988.  DIFFERENT GRAPHICS & ALPHA I/O (ALAN)
C     UPDATED         --FEBRUARY  1988.  UPDATED CYBER CONSTANTS (ALAN)
C     UPDATED         --JUNE      1989.  IBM-PC OS/2 & COMPAQ 386 CONSTANTS
C     UPDATED         --JUNE      1989.  INTEGER*2 (COMPAQ ERROR MESSAGE)
C     UPDATED         --JUNE      1990.  CODED COMMENTS FOR EASY "GLOBAL" EDIT
C                                        MOVE DATA AFTER EXECUTABLE
C     UPDATED         --AUGUST    1990.  (CONVEX, 4 DIFFERENT MODES, FROM CMLIB)
C     UPDATED         --APRIL     1992.  SAVE STATEMENTS
C     UPDATED         --APRIL     1992.  IHMOD1='386 '
C     UPDATED         --APRIL     1992.  ICOMPI='OTG '
C     UPDATED         --MAY       1992.  D.P. OVERFLOW PROBLEMS
C     UPDATED         --OCTOBER   1994.  FIX IBM-PC CONSTANTS
C     UPDATED         --OCTOBER   1996.  MS-FORTAN SUPPORT
C     UPDATED         --FEBRUARY  2006.  DOUBLE PRECISION VERSION
C                                        OF INTEL WINDOWS COMPILER
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IDEBUG
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION D2MACH(5)
C
CCCCC THE FOLLOWING 12 LINES WERE ADDED  APRIL 1992
C
      SAVE R2MACH
      SAVE ISMALL
      SAVE ILARGE
      SAVE IRIGHT
      SAVE IDIVER
      SAVE ILOG10
C
      SAVE D2MACH
      SAVE JSMALL
      SAVE JLARGE
      SAVE JRIGHT
      SAVE JDIVER
      SAVE JLOG10
C
      DIMENSION ISMALL(2)
      DIMENSION ILARGE(2)
      DIMENSION IRIGHT(2)
      DIMENSION IDIVER(2)
      DIMENSION ILOG10(2)
C
      DIMENSION JSMALL(4)
      DIMENSION JLARGE(4)
      DIMENSION JRIGHT(4)
      DIMENSION JDIVER(4)
      DIMENSION JLOG10(4)
C
      DIMENSION I2MACH(16)
C
      DIMENSION R2MACH(5)
C
      EQUIVALENCE (R2MACH(1),ISMALL(1))
      EQUIVALENCE (R2MACH(2),ILARGE(1))
      EQUIVALENCE (R2MACH(3),IRIGHT(1))
      EQUIVALENCE (R2MACH(4),IDIVER(1))
      EQUIVALENCE (R2MACH(5),ILOG10(1))
C
      EQUIVALENCE (D2MACH(1),JSMALL(1))
      EQUIVALENCE (D2MACH(2),JLARGE(1))
      EQUIVALENCE (D2MACH(3),JRIGHT(1))
      EQUIVALENCE (D2MACH(4),JDIVER(1))
      EQUIVALENCE (D2MACH(5),JLOG10(1))
C
C-----COMMON----------------------------------------------------------
C
C     THIS IS COMMON FOR MACHINE CONSTANTS
C
C               **************************
C               **************************
C               **  DEFINE THE INTEGER  **
C               **  MACHINE CONSTANTS   **
C               **  ARRAY               **
C               **************************
C               **************************
C
      COMMON /IMACH/I1MACH(16)
C
C               ********************************************************
C               ********************************************************
C               **  DEFINE THE REAL (SINGLE PRECISION FLOATING POINT)  *
C               **  MACHINE CONSTANTS                                  *
C               **  ARRAY                                              *
C               ********************************************************
C               ********************************************************
C
      COMMON /RMACH/R1MACH(5)
C
C               ***********************************
C               ***********************************
C               **  DEFINE THE DOUBLE PRECISION  **
C               **  MACHINE CONSTANTS            **
C               **  ARRAY                        **
C               ***********************************
C               ***********************************
C
      DOUBLE PRECISION D1MACH
C
      COMMON /DMACH/D1MACH(5)
C
C-----COMMON VARIABLES (GENERAL)--------------------------------------
C
C---------------------------------------------------------------------
C -----DATA STATEMENTS---------------------------------------------
CC
CC
CC    CONSTANTS FOR WINDOWS USING INTEL COMPILER
CC
#ifdef WIN32
C
      DATA I2MACH( 1) /    5 /
      DATA I2MACH( 2) /    6 /
      DATA I2MACH( 3) /    6 /
      DATA I2MACH( 4) /    0 /
      DATA I2MACH( 5) /   32 /
      DATA I2MACH( 6) /    4 /
      DATA I2MACH( 7) /    2 /
      DATA I2MACH( 8) /   31 /
      DATA I2MACH( 9) / 2147483647 /
      DATA I2MACH(10) /    2 /
      DATA I2MACH(11) /   24 /
      DATA I2MACH(12) / -125 /
      DATA I2MACH(13) /  127 /
      DATA I2MACH(14) /   53 /
      DATA I2MACH(15) / -1021 /
      DATA I2MACH(16) /  1023 /
C
      DATA R2MACH(3) / 5.9604645E-08 /
      DATA R2MACH(5) / 3.0102999E-01 /
      DATA D2MACH(3) / 1.1102230246252E-16  /
      DATA D2MACH(5) / 0.30102999563981195D0 /
C
      R2MACH(1)=TINY(1.0)
      R2MACH(2)=HUGE(1.0)
      R2MACH(4)=EPSILON(1.0)
      D2MACH(1)=TINY(1.0D0)
      D2MACH(2)=HUGE(1.0D0)
      D2MACH(4)=EPSILON(1.0D0)
#endif
CC
#ifdef LINUX
CC  FOLLOWING SECTION FOR LINUX.  SUPPORT CONDITIONAL COMPILATION.
CC
CC  FOR INTEGERS - SPECIFY "-DI32"  OR  "-DI64"  BASED ON WHETHER
CC  YOU ARE RUNNING ON A 32-BIT OR A 64-BIT MACHINE.
CC
      DATA I2MACH( 1) /    5 /
      DATA I2MACH( 2) /    6 /
      DATA I2MACH( 3) /    6 /
      DATA I2MACH( 4) /    0 /
      DATA I2MACH( 5) /   32 /
      DATA I2MACH( 6) /    4 /
      DATA I2MACH( 7) /    2 /
      DATA I2MACH( 8) /   31 /
      DATA I2MACH( 9) / 2147483647 /
      DATA I2MACH(10) /    2 /
      DATA I2MACH(11) /   53 /
      DATA I2MACH(12) / -1021 /
      DATA I2MACH(13) /  1024 /
      DATA I2MACH(14) /   113 /
      DATA I2MACH(15) / -16382 /
      DATA I2MACH(16) /  16384 /
C
      DATA R2MACH(1) / 1.175495E-38 /
      DATA R2MACH(2) / 3.402823E38 /
      DATA R2MACH(3) / 1.1920928955078E-7 /
      DATA R2MACH(4) / 2.3841857910156E-7 /
      DATA R2MACH(5) / 0.3010300 /
      DATA D2MACH(1) / 2.2250738585072D-307 /
      DATA D2MACH(2) / 1.7976931348623D308  /
      DATA D2MACH(3) / 1.1102230246252D-16  /
      DATA D2MACH(4) / 2.2204460492503D-16  /
      DATA D2MACH(5) / 0.30102999566398D0   /
#endif
C
C
C-----START POINT-----------------------------------------------------
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,51)
   51   FORMAT(1X)
        WRITE(IPR,52)
   52   FORMAT('***** AT THE BEGINNING OF INITMC--')
      ENDIF
C
C               *************************
C               **  COPY OVER INTEGER  **
C               **  MACHINE CONSTANTS  **
C               *************************
C
      DO100I=1,16
        I1MACH(I)=I2MACH(I)
  100 CONTINUE
C
C               ********************************************************
C               **  COPY OVER REAL (SINGLE PRECISION FLOATING POINT)  **
C               **  MACHINE CONSTANTS                                 **
C               ********************************************************
C
      DO200I=1,5
        R1MACH(I)=R2MACH(I)
  200 CONTINUE
C
CC              **********************************
CC              **  COPY OVER DOUBLE PRECISION  **
CC              **  MACHINE CONSTANTS           **
CC              **********************************
C
      DO300I=1,5
        D1MACH(I)=D2MACH(I)
  300 CONTINUE
C
      RETURN
      END
      SUBROUTINE MEAN(X,N,IPR,XMEAN,IDEBUG,IERROR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE SAMPLE MEAN
C              OF THE DATA IN THE INPUT VECTOR X.
C              THE SAMPLE MEAN = (SUM OF THE OBSERVATIONS)/N.
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 ARGUMENTS--XMEAN  = THE SINGLE PRECISION VALUE OF THE
C                                COMPUTED SAMPLE MEAN.
C     OUTPUT--THE COMPUTED SINGLE PRECISION VALUE OF THE SAMPLE MEAN.
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--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     REFERENCES--KENDALL AND STUART, THE ADVANCED THEORY OF
C                 STATISTICS, VOLUME 2, EDITION 1, 1961, PAGE 4.
C               --MOOD AND GRABLE, INTRODUCTION TO THE THEORY
C                 OF STATISTICS, EDITION 2, 1963, PAGE 146.
C               --DIXON AND MASSEY, INTRODUCTION TO STATISTICAL
C                 ANALYSIS, EDITION 2, 1957, PAGE 14.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY.
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--JUNE      1972.
C     UPDATED         --SEPTEMBER 1975.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --JUNE      1979.
C     UPDATED         --JULY      1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DN
      DOUBLE PRECISION DSUM
C
      DIMENSION X(*)
C
C-----START POINT-----------------------------------------------------
C
      ISUBN1='MEAN'
      ISUBN2='    '
      IERROR='NO'
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF MEAN--')
        WRITE(IPR,52)IDEBUG,N
   52   FORMAT('IDEBUG,N = ',A4,2X,I8)
        DO55I=1,N
          WRITE(IPR,56)I,X(I)
   56     FORMAT('I,X(I) = ',I8,G15.7)
   55   CONTINUE
      ENDIF
C
C               ********************
C               **  COMPUTE MEAN  **
C               ********************
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      AN=N
C
      IF(N.LT.1)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,111)
  111   FORMAT('***** ERROR IN MEAN--')
        WRITE(IPR,112)
  112   FORMAT('      THE NUMBER OF OBSERVATIONS IN THE RESPONSE ',
     1         'VARIABLE IS LESS THAN 1.')
        WRITE(IPR,117)N
  117   FORMAT('      THE NUMBER OF OBSERVATIONS = ',I8,'.')
        GOTO9000
      ELSEIF(N.EQ.1)THEN
        XMEAN=X(1)
        GOTO9000
      ENDIF
C
      HOLD=X(1)
      DO135I=2,N
      IF(X(I).NE.HOLD)GOTO139
  135 CONTINUE
      XMEAN=HOLD
      GOTO9000
  139 CONTINUE
C
C               *************************
C               **  STEP 2--           **
C               **  COMPUTE THE MEAN.  **
C               *************************
C
      DN=N
      DSUM=0.0D0
      DO200I=1,N
        DSUM=DSUM+DBLE(X(I))
  200 CONTINUE
      XMEAN=REAL(DSUM/DN)
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT.      **
C               *****************
C
 9000 CONTINUE
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF MEAN--')
        WRITE(IPR,9012)IERROR,XMEAN
 9012   FORMAT('IERROR,XMEAN = ',A4,2X,G15.7)
      ENDIF
C
      RETURN
      END
      SUBROUTINE NORPPF(P,PPF,IPR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE FOR THE NORMAL (GAUSSIAN)
C              DISTRIBUTION WITH MEAN = 0 AND STANDARD DEVIATION = 1.
C              THIS DISTRIBUTION IS DEFINED FOR ALL X AND HAS
C              THE PROBABILITY DENSITY FUNCTION
C              F(X) = (1/SQRT(2*PI))*EXP(-X*X/2).
C              NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION
C              IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE
C              DISTRIBUTION FUNCTION OF THE DISTRIBUTION.
C     INPUT  ARGUMENTS--P      = THE SINGLE PRECISION VALUE
C                                (BETWEEN 0.0 AND 1.0)
C                                AT WHICH THE PERCENT POINT
C                                FUNCTION IS TO BE EVALUATED.
C     OUTPUT ARGUMENTS--PPF    = THE SINGLE PRECISION PERCENT
C                                POINT FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION PERCENT POINT
C             FUNCTION VALUE PPF.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
C     RESTRICTIONS--P SHOULD BE BETWEEN 0.0 AND 1.0, EXCLUSIVELY.
C     OTHER DATAPAC   SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, LOG.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     REFERENCES--ODEH AND EVANS, THE PERCENTAGE POINTS
C                 OF THE NORMAL DISTRIBUTION, ALGORTIHM 70,
C                 APPLIED STATISTICS, 1974, PAGES 96-97.
C               --EVANS, ALGORITHMS FOR MINIMAL DEGREE
C                 POLYNOMIAL AND RATIONAL APPROXIMATION,
C                 M. SC. THESIS, 1972, UNIVERSITY
C                 OF VICTORIA, B. C., CANADA.
C               --HASTINGS, APPROXIMATIONS FOR DIGITAL
C                 COMPUTERS, 1955, PAGES 113, 191, 192.
C               --NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 933, FORMULA 26.2.23.
C               --FILLIBEN, SIMPLE AND ROBUST LINEAR ESTIMATION
C                 OF THE LOCATION PARAMETER OF A SYMMETRIC
C                 DISTRIBUTION (UNPUBLISHED PH.D. DISSERTATION,
C                 PRINCETON UNIVERSITY), 1969, PAGES 21-44, 229-231.
C               --FILLIBEN, 'THE PERCENT POINT FUNCTION',
C                 (UNPUBLISHED MANUSCRIPT), 1970, PAGES 28-31.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--1, 1970, PAGES 40-111.
C               --THE KELLEY STATISTICAL TABLES, 1948.
C               --OWEN, HANDBOOK OF STATISTICAL TABLES,
C                 1962, PAGES 3-16.
C               --PEARSON AND HARTLEY, BIOMETRIKA TABLES
C                 FOR STATISTICIANS, VOLUME 1, 1954,
C                 PAGES 104-113.
C     COMMENTS--THE CODING AS PRESENTED BELOW
C               IS ESSENTIALLY IDENTICAL TO THAT
C               PRESENTED BY ODEH AND EVANS
C               AS ALGORTIHM 70 OF APPLIED STATISTICS.
C               THE PRESENT AUTHOR HAS MODIFIED THE
C               ORIGINAL ODEH AND EVANS CODE WITH ONLY
C               MINOR STYLISTIC CHANGES.
C             --AS POINTED OUT BY ODEH AND EVANS
C               IN APPLIED STATISTICS,
C               THEIR ALGORITHM REPRESENTES A
C               SUBSTANTIAL IMPROVEMENT OVER THE
C               PREVIOUSLY EMPLOYED
C               HASTINGS APPROXIMATION FOR THE
C               NORMAL PERCENT POINT FUNCTION--
C               THE ACCURACY OF APPROXIMATION
C               BEING IMPROVED FROM 4.5*(10**-4)
C               TO 1.5*(10**-8).
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--JUNE      1972.
C     UPDATED         --SEPTEMBER 1975.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --OCTOBER   1976.
C     UPDATED         --MAY       1982.
C
C-----DATA STATEMENTS-------------------------------------------------
C
      DATA P0,P1,P2,P3,P4
     1/-.322232431088,-1.0,
     1 -.342242088547,-.204231210245E-1,
     1 -.453642210148E-4/
      DATA Q0,Q1,Q2,Q3,Q4
     1/.993484626060E-1,.588581570495,
     1 .531103462366,.103537752850,
     1 .38560700634E-2/
C
C-----START POINT-----------------------------------------------------
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      IF(P.LE.0.0 .OR. P.GE.1.0)THEN
        WRITE(IPR,1)
    1   FORMAT('***** ERROR IN NORPPF--')
        WRITE(IPR,2)
    2   FORMAT('      THE FIRST ARGUMENT TO NORPPF IS OUTSIDE ',
     1         'THE ALLOWABLE (0,1) INTERVAL.')
        WRITE(IPR,46)P
   46   FORMAT('      THE VALUE OF THE ARGUMENT IS ',G15.7,'.')
        GOTO9000
      ENDIF
C
C               ************************************************
C               **  STEP 2--                                  **
C               **  COMPUTE THE PERCENT POINT FUNCTION VALUE  **
C               ************************************************
C
      IF(P.EQ.0.5)THEN
        PPF=0.0
        GOTO9000
      ENDIF
C
      PHOLD=P
      R=P
      IF(PHOLD.GT.0.5)R=1.0-R
      T=SQRT(-2.0*LOG(R))
      ANUM=((((T*P4+P3)*T+P2)*T+P1)*T+P0)
      ADEN=((((T*Q4+Q3)*T+Q2)*T+Q1)*T+Q0)
      PPF=T+(ANUM/ADEN)
      IF(PHOLD.LT.0.5)PPF=-PPF
C
 9000 CONTINUE
      RETURN
      END
      SUBROUTINE PERCEN(P100,X,N,IPR,IQUAME,XTEMP,MAXNXT,
     1                  XPERC,IDEBUG,IERROR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE SAMPLE PERCENTILE
C              OF THE DATA IN THE INPUT VECTOR X.
C     INPUT  ARGUMENTS--P100   = THE SINGLE PRECISION PERCENTAGE
C                                (BETWEEN 0 AND 100)
C                     --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 ARGUMENTS--XPERC  = THE SINGLE PRECISION VALUE OF THE
C                                COMPUTED SAMPLE PERCENTILE.
C     OUTPUT--THE COMPUTED SINGLE PRECISION VALUE OF THE
C             SAMPLE PERCENTILE.
C     OTHER DATAPAC   SUBROUTINES NEEDED--SORT.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     REFERENCES--HYNDMAN AND FAN (NOVEMBER 1996), "SAMPLE QUANTILES IN
C                 STATISTICAL PACKAGES", THE AMERICAN STATISTICIAN,
C                 VOL. 50, NO. 4, PP. 361-365.
C     WRITTEN BY--JAMES J. FILLIBEN/ALAN HECKERT
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     VERSION NUMBER--87.11
C     ORIGINAL VERSION--SEPTEMBER 1987.
C     UPDATED         --FEBRUARY  2015. SUPPORT "R8" AND "R7" METHODS
C     UPDATED         --DECEMBER  2017. SUPPORT "AVERAGE" METHOD
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IQUAME
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
      DIMENSION XTEMP(*)
C
C---------------------------------------------------------------------
C
C-----START POINT-----------------------------------------------------
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF PERCEN--')
        WRITE(IPR,53)P100,N,IDEBUG
   53   FORMAT('P100,N,IDEBUG, = ',G15.7,I8,2X,A4)
        DO55I=1,N
          WRITE(IPR,56)I,X(I)
   56     FORMAT('I,X(I) = ',I8,G15.7)
   55   CONTINUE
      ENDIF
C
      ISUBN1='PERC'
      ISUBN2='EN  '
      IERROR='NO'
C
      NI=0
      NIP1=0
      ANI=0.0
      A2NI=0.0
      REM=0.0
C
C               ******************************
C               **  COMPUTE PERCENTILE      **
C               ******************************
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      AN=N
C
      IF(P100.LT.0.0 .OR. P100.GT.100.0)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,111)
  111   FORMAT('***** ERROR IN PERCENTILE--')
        WRITE(IPR,112)
  112   FORMAT('      THE INPUT PERCENTAGE FOR WHICH THE PERCENTILE ',
     1         'IS TO BE')
        WRITE(IPR,115)
  115   FORMAT('      COMPUTED MUST BE LARGER THAN 0 AND SMALLER ',
     1         'THAN 100.')
        WRITE(IPR,117)P100
  117   FORMAT('      SUCH WAS NOT THE CASE HERE.  THE INPUT ',
     1         'PERCENTAGE = ',G15.7)
        GOTO9000
      ENDIF
C
      IF(N.LT.1 .OR. N.GT.MAXNXT)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,111)
        WRITE(IPR,122)
  122   FORMAT('      THE INPUT NUMBER OF OBSERVATIONS FOR THE ',
     1         'RESPONSE VARIABLE')
        WRITE(IPR,125)MAXNXT
  125   FORMAT('      MUST BE BETWEEN 1 AND ',I8,' (INCLUSIVELY).')
        WRITE(IPR,127)N
  127   FORMAT('      THE INPUT NUMBER OF OBSERVATIONS = ',I8,'.')
        GOTO9000
      ENDIF
C
      IF(N.EQ.1)THEN
        XPERC=X(1)
        GOTO9000
      ENDIF
C
C     CASE WHERE ALL RESPONSE VALUES ARE THE SAME
C
      HOLD=X(1)
      DO145I=2,N
        IF(X(I).NE.HOLD)GOTO149
  145 CONTINUE
      XPERC=HOLD
      GOTO9000
  149 CONTINUE
C
C               ***********************************
C               **  STEP 2--                     **
C               **  COMPUTE THE PERCENTILE.      **
C               ***********************************
C
      CALL SORT(X,N,XTEMP,IPR)
C
      P=P100/100.0
C
C     NOTE 2015/02: THE HYNDMAN AND FAN ARTICLE DISCUSSES EIGHT
C                   DEFINITIONS OF A PERCENTILE IN RELATION TO SIX
C                   DESIRABLE PROPERTIES.  THE DEFAULT DATAPLOT
C                   ALGORITHM USES THE "R6" DEFINITION.  THEIR
C                   RECOMMENDATION FOR A "STANDARD" DEFINITION IS
C                   "R8".  ALLOW THE USER TO REQUEST THE R8
C                   DEFINITION (SPECIFIED IN IQUAME VARIABLE).
C                   THE R7 METHOD IS ALSO SUPPORTED (THIS IS THE
C                   DEFAULT METHOD FOR R AND EXCEL).
C
C                   ALSO, FOR DEFAULT DEFINITION, LOOK FOR MINIMUM
C                   AND MAXIMUM VALUES OF P100 THAT ARE SET TO THE
C                   DATA MINIMUM AND MAXIMUM, RESPECTIVELY.
C
      IF(IQUAME.EQ.'R8')THEN
        AVAL1=2.0/3.0
        AVAL2=1.0/3.0
        ALOW=AVAL1/(AN+AVAL2)
        AHIGH=(AN-AVAL2)/(AN+AVAL2)
        IF(P.LT.ALOW)THEN
          XPERC=XTEMP(1)
        ELSEIF(P.GE.AHIGH)THEN
          XPERC=XTEMP(N)
        ELSE
          ANI=P*(AN+AVAL2) + AVAL2
          NI=ANI
          A2NI=NI
          REM=ANI-A2NI
          NIP1=NI+1
          IF(NI.LE.1)NI=1
          IF(NI.GE.N)NI=N
          IF(NIP1.LE.1)NIP1=1
          IF(NIP1.GE.N)NIP1=N
          XPERC=XTEMP(NI) + REM*(XTEMP(NIP1) - XTEMP(NI))
        ENDIF
      ELSEIF(IQUAME.EQ.'R7')THEN
        AHIGH=1.0
        IF(P.GE.AHIGH)THEN
          XPERC=XTEMP(N)
        ELSE
          AVAL1=1.0
          ANI=P*(AN-AVAL1) + AVAL1
          NI=ANI
          A2NI=NI
          REM=ANI-A2NI
          NIP1=NI+1
          IF(NI.LE.1)NI=1
          IF(NI.GE.N)NI=N
          IF(NIP1.LE.1)NIP1=1
          IF(NIP1.GE.N)NIP1=N
          XPERC=XTEMP(NI) + REM*(XTEMP(NIP1) - XTEMP(NI))
        ENDIF
      ELSEIF(IQUAME.EQ.'AVER')THEN
C
C        AVERAGE METHOD - TAKE AVERAGE OF TWO CLOSEST POINTS, NOT
C        RECOMMENDED, BUT THIS APPEARS TO BE METHOD USED BY BONETT
C        FOR CONFIDENCE INTERVALS FOR COEFFCIENT OF QUARTILE
C        DISPERSION.  INCLUDE SO THAT WE CAN MATCH HIS RESULTS.
C
        AINDX=P*AN
        IINDX=INT(AINDX)
        IF(IINDX.EQ.N)THEN
          XPERC=XTEMP(IINDX)
        ELSEIF(IINDX.LT.1)THEN
          XPERC=XTEMP(1)
        ELSE
          XPERC=(XTEMP(IINDX) + XTEMP(IINDX+1))/2.0
        ENDIF
      ELSE
        ALOW=1.0/(AN+1.0)
        AHIGH=AN/(AN+1.0)
        IF(P.LT.ALOW)THEN
          XPERC=XTEMP(1)
        ELSEIF(P.GE.AHIGH)THEN
          XPERC=XTEMP(N)
        ELSE
          ANI=P*(AN+1.0)
          NI=ANI
          A2NI=NI
          REM=ANI-A2NI
          NIP1=NI+1
          IF(NI.LE.1)NI=1
          IF(NI.GE.N)NI=N
          IF(NIP1.LE.1)NIP1=1
          IF(NIP1.GE.N)NIP1=N
CCCCC     BUG FIX.  WEIGHTS IN WRONG ORDER!  NOVEMBER 1998.
CCCCC     XPERC=REM*XTEMP(NI)+(1.0-REM)*XTEMP(NIP1)
          XPERC=(1.0-REM)*XTEMP(NI)+REM*XTEMP(NIP1)
        ENDIF
      ENDIF
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT.      **
C               *****************
C
 9000 CONTINUE
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF PERCEN--')
        WRITE(IPR,9013)IERROR,N,P100,P,XPERC
 9013   FORMAT('IERROR,N,P100,P,XPERC = ',A4,2X,I8,3G15.7)
        WRITE(IPR,9014)NI,NIP1,ANI,A2NI,REM
 9014   FORMAT('NI,NIP1,ANI,A2NI,REM = ',2I8,3G15.7)
      ENDIF
C
      RETURN
      END
      SUBROUTINE SD(X,N,IPR,XSD,IDEBUG,IERROR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE SAMPLE STANDARD DEVIATION
C              (WITH DENOMINATOR N-1) OF THE DATA IN THE INPUT VECTOR X.
C              THE SAMPLE STANDARD DEVIATION = SQRT((THE SUM OF THE
C              SQUARED DEVIATIONS ABOUT THE SAMPLE MEAN)/(N-1)).
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 ARGUMENTS--XSD    = THE SINGLE PRECISION VALUE OF THE
C                                COMPUTED SAMPLE STANDARD DEVIATION.
C     OUTPUT--THE COMPUTED SINGLE PRECISION VALUE OF THE
C             SAMPLE STANDARD DEVIATION (WITH DENOMINATOR N-1).
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--DSQRT.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     REFERENCES--SNEDECOR AND COCHRAN, STATISTICAL METHODS,
C                 EDITION 6, 1967, PAGE 44.
C               --DIXON AND MASSEY, INTRODUCTION TO STATISTICAL
C                 ANALYSIS, EDITION 2, 1957, PAGES 19, 76.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATION INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATION INSTITUTE OF STANDARDS AND TECHNOLOGY.
C     LANGUAGE--ANSI FORTRAN (1966)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--JUNE      1972.
C     UPDATED         --SEPTEMBER 1975.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --AUGUST    1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DN
      DOUBLE PRECISION DX
      DOUBLE PRECISION DSUM
      DOUBLE PRECISION DMEAN
      DOUBLE PRECISION DVAR
      DOUBLE PRECISION DSD
C
      DIMENSION X(*)
C
C-----START POINT-----------------------------------------------------
C
      ISUBN1='SD  '
      ISUBN2='    '
      IERROR='NO'
C
      XSD=0.0
      DMEAN=0.0D0
      AN=N
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF SD--')
        WRITE(IPR,52)IDEBUG,N
   52   FORMAT('IDEBUG,N = ',A4,2X,I8)
        DO55I=1,N
          WRITE(IPR,56)I,X(I)
   56     FORMAT('I,X(I) = ',I8,G15.7)
   55   CONTINUE
      ENDIF
C
C               **********************************
C               **  COMPUTE STANDARD DEVIATION  **
C               **********************************
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      IF(N.LT.1)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,111)
  111   FORMAT('***** ERROR IN SD--')
        WRITE(IPR,112)
  112   FORMAT('      THE NUMBER OF OBSERVATIONS FOR THE RESPONSE')
        WRITE(IPR,115)
  115   FORMAT('      VARIABLE IS LESS THAN 1.')
        WRITE(IPR,117)N
  117   FORMAT('      THE INPUT NUMBER OF OBSERVATIONS = ',I8,'.')
        GOTO9000
      ELSEIF(N.EQ.1)THEN
        GOTO9000
      ENDIF
C
      HOLD=X(1)
      DO135I=2,N
        IF(X(I).NE.HOLD)GOTO139
  135 CONTINUE
      GOTO9000
  139 CONTINUE
C
C               ***************************************
C               **  STEP 2--                         **
C               **  COMPUTE THE STANDARD DEVIATION.  **
C               ***************************************
C
      DN=N
      DSUM=0.0D0
      DO200I=1,N
        DX=X(I)
        DSUM=DSUM+DX
  200 CONTINUE
      DMEAN=DSUM/DN
C
      DSUM=0.0D0
      DO300I=1,N
        DX=X(I)
        DSUM=DSUM+(DX-DMEAN)**2
  300 CONTINUE
      DVAR=DSUM/(DN-1.0D0)
      DSD=0.0D0
      IF(DVAR.GT.0.0D0)DSD=DSQRT(DVAR)
      XSD=DSD
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT.      **
C               *****************
C
 9000 CONTINUE
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF SD--')
        WRITE(IPR,9012)IERROR,DMEAN,XSD
 9012   FORMAT('IERROR,DMEAN,XSD = ',A4,2X,2G15.7)
      ENDIF
C
      RETURN
      END
      SUBROUTINE SORT(X,N,Y,IPR)
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 (1977)
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  QUICKSORT.
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         QUICKSORT       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 ( QUICKSORT 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 DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATION INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--JUNE      1972.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --JUNE      1981.
C     UPDATED         --AUGUST    1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
      DIMENSION Y(*)
C
      DIMENSION IU(36)
      DIMENSION IL(36)
C
C-----START POINT-----------------------------------------------------
C
      ISUBN1='SORT'
      ISUBN2='    '
      IERROR='NO'
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF SORT--')
        WRITE(IPR,52)IDEBUG,N
   52   FORMAT('IDEBUG,N = ',A4,2X,I8)
        DO55I=1,N
          WRITE(IPR,56)I,X(I)
   56     FORMAT('I,X(I) = ',I8,G15.7)
   55   CONTINUE
      ENDIF
C
C               ************************
C               **  SORT THE VALUES.  **
C               ************************
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      IF(N.LT.1)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,111)
  111   FORMAT('***** ERROR IN SORT--')
        WRITE(IPR,112)
  112   FORMAT('      THE NUMBER OF OBSERVATIONS TO BE SORTED IS ',
     1         'LESS THAN 1.')
        WRITE(IPR,118)N
  118   FORMAT('      THE NUMBER OF OBSERVATIONS IS ',I8)
        GOTO9000
      ELSEIF(N.EQ.1)THEN
        Y(1)=X(1)
        GOTO9000
      ENDIF
C
C               *******************************************
C               **  STEP 2--                             **
C               **  COPY THE VECTOR X INTO THE VECTOR Y  **
C               *******************************************
C
      DO137I=1,N
        Y(I)=X(I)
  137 CONTINUE
C
C     CHECK FOR CASE WHERE ALL DATA HAVE THE SAME VALUE (SO
C     ALREADY SORTED)
C
      HOLD=X(1)
      DO135I=2,N
        IF(X(I).NE.HOLD)GOTO139
  135 CONTINUE
      GOTO9000
C
  139 CONTINUE
C
C               **********************************************************
C               **  STEP 3--                                            **
C               **  CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED  **
C               **********************************************************
C
      NM1=N-1
      DO250I=1,NM1
        IP1=I+1
        IF(Y(I).LE.Y(IP1))GOTO250
        GOTO290
  250 CONTINUE
      GOTO9000
  290 CONTINUE
C
C               ***************************
C               **  STEP 4--             **
C               **  CARRY OUT THE SORT.  **
C               ***************************
C
      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)GOTO9000
      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
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT.      **
C               *****************
C
 9000 CONTINUE
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF SORT--')
        WRITE(IPR,9012)IERROR
 9012   FORMAT('IERROR = ',A4)
        DO9015I=1,N
          WRITE(IPR,9016)I,X(I),Y(I)
 9016     FORMAT('I,X(I),Y(I) = ',I8,2G15.7)
 9015   CONTINUE
      ENDIF
C
      RETURN
      END
      SUBROUTINE SORTC(X,Y,N,XS,YC,IPR)
C
C     PURPOSE--THIS SUBROUTINE SORTS (IN ASCENDING ORDER)
C              THE N ELEMENTS OF THE SINGLE PRECISION VECTOR X,
C              PUTS THE RESULTING N SORTED VALUES INTO THE
C              SINGLE PRECISION VECTOR XS,
C              REARRANGES THE ELEMENTS OF THE VECTOR Y
C              (ACCORDING TO THE SORT ON X),
C              AND PUTS THE REARRANGED Y VALUES
C              INTO THE SINGLE PRECISION VECTOR YC.
C              THIS SUBROUTINE GIVES THE DATA ANALYST
C              THE ABILITY TO SORT ONE DATA VECTOR
C              WHILE 'CARRYING ALONG' THE ELEMENTS
C              OF A SECOND DATA VECTOR.
C     INPUT  ARGUMENTS--X      = THE SINGLE PRECISION VECTOR OF
C                                OBSERVATIONS TO BE SORTED.
C                     --Y      = THE SINGLE PRECISION VECTOR OF
C                                OBSERVATIONS TO BE 'CARRIED ALONG',
C                                THAT IS, TO BE REARRANGED ACCORDING
C                                TO THE SORT ON X.
C                     --N      = THE INTEGER NUMBER OF OBSERVATIONS
C                                IN THE VECTOR X.
C     OUTPUT ARGUMENTS--XS     = THE SINGLE PRECISION VECTOR
C                                INTO WHICH THE SORTED DATA VALUES
C                                FROM X WILL BE PLACED.
C                     --YC     = THE SINGLE PRECISION VECTOR
C                                INTO WHICH THE REARRANGED
C                                (ACCORDING TO THE SORT OF THE
C                                VECTOR X) VALUES OF THE VECTOR Y
C                                WILL BE PLACED.
C     OUTPUT--THE SINGLE PRECISION VECTOR XS
C             CONTAINING THE SORTED
C             (IN ASCENDING ORDER) VALUES
C             OF THE SINGLE PRECISION VECTOR X, AND
C             THE SINGLE PRECISION VECTOR YC
C             CONTAINING THE REARRANGED
C             (ACCORDING TO THE SORT ON X)
C             VALUES OF THE VECTOR Y.
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 (1977)
C     COMMENT--THE SMALLEST ELEMENT OF THE VECTOR X
C              WILL BE PLACED IN THE FIRST POSITION
C              OF THE VECTOR XS,
C              THE SECOND SMALLEST ELEMENT IN THE VECTOR X
C              WILL BE PLACED IN THE SECOND POSITION
C              OF THE VECTOR XS,
C              ETC.
C     COMMENT--THE ELEMENT IN THE VECTOR Y CORRESPONDING
C              TO THE SMALLEST ELEMENT IN X
C              WILL BE PLACED IN THE FIRST POSITION
C              OF THE VECTOR YC,
C              THE ELEMENT IN THE VECTOR Y CORRESPONDING
C              TO THE SECOND SMALLEST ELEMENT IN X
C              WILL BE PLACED IN THE SECOND POSITION
C              OF THE VECTOR YC,
C              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 SORTC(X,Y,N,X,YC)
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 DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATION INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--JUNE      1972.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --JUNE      1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DIMENSION X(*),Y(*),XS(*),YC(*)
      DIMENSION IU(36),IL(36)
C
C-----START POINT-----------------------------------------------------
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)THEN
        WRITE(IPR,15)
   15   FORMAT('***** ERROR IN SORTC--')
        WRITE(IPR,16)
   16   FORMAT('      THE NUMBER OF OBSERVATIONS IS NON-POSITIVE.')
        WRITE(IPR,47)N
   47   FORMAT('      THE NUMBER OF OBSERVATIONS IS ',I8)
        GOTO9000
      ELSEIF(N.EQ.1)THEN
        XS(1)=X(1)
        YC(1)=Y(1)
        GOTO9000
      ENDIF
C
C     COPY THE VECTOR X INTO THE VECTOR XS AND THE VECTOR Y INTO YC
C
      DO100I=1,N
        XS(I)=X(I)
        YC(I)=Y(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(XS(I).LE.XS(IP1))GOTO200
      GOTO250
  200 CONTINUE
      GOTO9000
  250 M=1
      I=1
      J=N
  305 IF(I.GE.J)GOTO370
  310 K=I
      MID=(I+J)/2
      AMED=XS(MID)
      BMED=YC(MID)
      IF(XS(I).LE.AMED)GOTO320
      XS(MID)=XS(I)
      YC(MID)=YC(I)
      XS(I)=AMED
      YC(I)=BMED
      AMED=XS(MID)
      BMED=YC(MID)
  320 L=J
      IF(XS(J).GE.AMED)GOTO340
      XS(MID)=XS(J)
      YC(MID)=YC(J)
      XS(J)=AMED
      YC(J)=BMED
      AMED=XS(MID)
      BMED=YC(MID)
      IF(XS(I).LE.AMED)GOTO340
      XS(MID)=XS(I)
      YC(MID)=YC(I)
      XS(I)=AMED
      YC(I)=BMED
      AMED=XS(MID)
      BMED=YC(MID)
      GOTO340
  330 XS(L)=XS(K)
      YC(L)=YC(K)
      XS(K)=TX
      YC(K)=TY
  340 L=L-1
      IF(XS(L).GT.AMED)GOTO340
      TX=XS(L)
      TY=YC(L)
  350 K=K+1
      IF(XS(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)GOTO9000
      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=XS(I+1)
      BMED=YC(I+1)
      IF(XS(I).LE.AMED)GOTO390
      K=I
  395 XS(K+1)=XS(K)
      YC(K+1)=YC(K)
      K=K-1
      IF(AMED.LT.XS(K))GOTO395
      XS(K+1)=AMED
      YC(K+1)=BMED
      GOTO390
C
 9000 CONTINUE
      RETURN
      END
      SUBROUTINE TPPF(P,ANU,PPF,IPR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT
C              FUNCTION VALUE FOR THE STUDENT'S T DISTRIBUTION
C              WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU.
C              THE STUDENT'S T DISTRIBUTION USED
C              HEREIN IS DEFINED FOR ALL X,
C              AND ITS PROBABILITY DENSITY FUNCTION IS GIVEN
C              IN THE REFERENCES BELOW.
C              NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION
C              IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE
C              DISTRIBUTION FUNCTION OF THE DISTRIBUTION.
C     INPUT  ARGUMENTS--P      = THE SINGLE PRECISION VALUE
C                                (BETWEEN 0.0 (EXCLUSIVELY)
C                                AND 1.0 (EXCLUSIVELY))
C                                AT WHICH THE PERCENT POINT
C                                FUNCTION IS TO BE EVALUATED.
C                     --NU     = THE INTEGER NUMBER OF DEGREES
C                                OF FREEDOM.
C                                NU SHOULD BE POSITIVE.
C     OUTPUT ARGUMENTS--PPF    = THE SINGLE PRECISION PERCENT
C                                POINT FUNCTION VALUE.
C     OUTPUT--THE SINGLE PRECISION PERCENT POINT FUNCTION .
C             VALUE PPF FOR THE STUDENT'S T DISTRIBUTION
C             WITH DEGREES OF FREEDOM PARAMETER = NU.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
C     RESTRICTIONS--NU SHOULD BE A POSITIVE INTEGER VARIABLE.
C                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--NORPPF.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--DSIN, DCOS, DSQRT, DATAN.
C     MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     COMMENT--FOR NU = 1 AND NU = 2, THE PERCENT POINT FUNCTION
C              FOR THE T DISTRIBUTION EXISTS IN SIMPLE CLOSED FORM
C              AND SO THE COMPUTED PERCENT POINTS ARE EXACT.
C            --FOR OTHER SMALL VALUES OF NU (NU BETWEEN 3 AND 6,
C              INCLUSIVELY), THE APPROXIMATION
C              OF THE T PERCENT POINT BY THE FORMULA
C              GIVEN IN THE REFERENCE BELOW IS AUGMENTED
C              BY 3 ITERATIONS OF NEWTON'S METHOD FOR
C              ROOT DETERMINATION.
C              THIS IMPROVES THE ACCURACY--ESPECIALLY FOR
C              VALUES OF P NEAR 0 OR 1.
C     REFERENCES--NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY APPLIED MATHMATICS
C                 SERIES 55, 1964, PAGE 949, FORMULA 26.7.5.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--2, 1970, PAGE 102,
C                 FORMULA 11.
C               --FEDERIGHI, 'EXTENDED TABLES OF THE
C                 PERCENTAGE POINTS OF STUDENT'S T
C                 DISTRIBUTION, JOURNAL OF THE
C                 AMERICAN STATISTICAL ASSOCIATION,
C                 1969, PAGES 683-688.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGES 120-123.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-921-3651
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--OCTOBER   1975.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --DECEMBER  1981.
C     UPDATED         --MAY       1982.
C     UPDATED         --OCTOBER   2006. SUPPORT FOR FRACTIONAL
C                                       DEGREES OF FREEDOM
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION PI
      DOUBLE PRECISION SQRT2
      DOUBLE PRECISION DP
      DOUBLE PRECISION DNU
      DOUBLE PRECISION TERM1,TERM2,TERM3,TERM4,TERM5
      DOUBLE PRECISION DPPFN
      DOUBLE PRECISION DPPF,DCON,DARG,Z,S,C
      DOUBLE PRECISION B21
      DOUBLE PRECISION B31,B32,B33,B34
      DOUBLE PRECISION B41,B42,B43,B44,B45
      DOUBLE PRECISION B51,B52,B53,B54,B55,B56
      DOUBLE PRECISION D1,D3,D5,D7,D9
C
      DOUBLE PRECISION DC
      DOUBLE PRECISION DALPHA
      DOUBLE PRECISION EPS
      DOUBLE PRECISION SIG
      DOUBLE PRECISION ZERO
      DOUBLE PRECISION DMEAN
      DOUBLE PRECISION DSD
      DOUBLE PRECISION XL
      DOUBLE PRECISION XR
      DOUBLE PRECISION XINC
      DOUBLE PRECISION X
      DOUBLE PRECISION FXL
      DOUBLE PRECISION FXR
      DOUBLE PRECISION P1
      DOUBLE PRECISION FCS
      DOUBLE PRECISION XRML
      DOUBLE PRECISION DCDF
      DOUBLE PRECISION CDFL
      DOUBLE PRECISION CDFR
      DOUBLE PRECISION DBETAI
C
      EXTERNAL DBETAI
C
C-----DATA STATEMENTS-------------------------------------------------
C
      DATA PI/3.14159265358979D0/
      DATA SQRT2/1.414213562D0/
      DATA B21/0.25D0/
      DATA B31,B32,B33,B34/0.01041666666667D0,5.0D0,16.0D0,3.0D0/
      DATA B41,B42,B43,B44,B45/0.00260416666667D0,3.0D0,19.0D0,17.0D0,
     1                         -15.0D0/
      DATA B51,B52,B53,B54,B55,B56/0.00001085069444D0,79.0D0,776.0D0,
     1                             1482.0D0,-1920.0D0,-945.0D0/
C
      DATA EPS /0.0001D0/
      DATA SIG /1.0D-7/
      DATA ZERO /0.0D0/
C
C-----START POINT-----------------------------------------------------
C
      S=0.0D0
      C=0.0D0
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(P.LE.0.0.OR.P.GE.1.0)THEN
        WRITE(IPR,1)
    1   FORMAT('***** ERROR IN TPPF--')
        WRITE(IPR,2)
    2   FORMAT('      THE FIRST ARGUMENT TO TPPF IS OUTSIDE THE ',
     1         'ALLOWABLE (0,1) INTERVAL')
        WRITE(IPR,46)P
   46   FORMAT('      THE VALUE OF THE ARGUMENT IS ',G15.7)
        GOTO9000
      ENDIF
C
      IF(ANU.LT.1.0)THEN
        WRITE(IPR,1)
        WRITE(IPR,11)
   11   FORMAT('      THE DEGREES OF FREEDOM ARGUMENT TO TPPF ',
     1         'IS LESS THAN 1')
        WRITE(IPR,46)ANU
        PPF=0.0
        GOTO9000
      ENDIF
C
      IF(P.EQ.0.5)THEN
        PPF=0.0
        GOTO9000
      ENDIF
C
      NU=INT(ANU)
      IF(ABS(ANU-REAL(NU)).GT.0.000001)GOTO8000
C
      DNU=NU
      DP=P
      MAXIT=5
C
      IF(NU.EQ.1)THEN
C
C       TREAT THE NU = 1 (CAUCHY) CASE
C
        DARG=PI*DP
        PPF=-DCOS(DARG)/DSIN(DARG)
        GOTO9000
      ELSEIF(NU.EQ.2)THEN
C
C       TREAT THE NU = 2 CASE
C
        TERM1=SQRT2/2.0D0
        TERM2=2.0D0*DP-1.0D0
        TERM3=DSQRT(DP*(1.0D0-DP))
        PPF=TERM1*TERM2/TERM3
        GOTO9000
      ELSE
C
C       TREAT THE NU GREATER THAN OR EQUAL TO 3 CASE
C
        CALL NORPPF(P,PPFN,IPR)
        DPPFN=PPFN
        D1=DPPFN
        D3=DPPFN**3
        D5=DPPFN**5
        D7=DPPFN**7
        D9=DPPFN**9
        TERM1=D1
        TERM2=B21*(D3+D1)/DNU
        TERM3=B31*(B32*D5+B33*D3+B34*D1)/(DNU**2)
        TERM4=B41*(B42*D7+B43*D5+B44*D3+B45*D1)/(DNU**3)
        TERM5=B51*(B52*D9+B53*D7+B54*D5+B55*D3+B56*D1)/(DNU**4)
        DPPF=TERM1+TERM2+TERM3+TERM4+TERM5
        PPF=DPPF
        IF(NU.GE.7)GOTO9000
        IF(NU.EQ.3)THEN
C
C         AUGMENT THE RESULTS FOR THE NU = 3 CASE
C
          DCON=PI*(DP-0.5D0)
          DARG=DPPF/DSQRT(DNU)
          Z=DATAN(DARG)
          DO350IPASS=1,MAXIT
            S=DSIN(Z)
            C=DCOS(Z)
            Z=Z-(Z+S*C-DCON)/(2.0D0*C*C)
  350     CONTINUE
          PPF=DSQRT(DNU)*S/C
        ELSEIF(NU.EQ.4)THEN
C
C         AUGMENT THE RESULTS FOR THE NU = 4 CASE
C
          DCON=2.0D0*(DP-0.5D0)
          DARG=DPPF/DSQRT(DNU)
          Z=DATAN(DARG)
          DO450IPASS=1,MAXIT
            S=DSIN(Z)
            C=DCOS(Z)
            Z=Z-((1.0D0+0.5D0*C*C)*S-DCON)/(1.5D0*C*C*C)
  450     CONTINUE
          PPF=DSQRT(DNU)*S/C
        ELSEIF(NU.EQ.5)THEN
C
C         AUGMENT THE RESULTS FOR THE NU = 5 CASE
C
          DCON=PI*(DP-0.5D0)
          DARG=DPPF/DSQRT(DNU)
          Z=DATAN(DARG)
          DO550IPASS=1,MAXIT
            S=DSIN(Z)
            C=DCOS(Z)
            Z=Z-(Z+(C+(2.0D0/3.0D0)*C*C*C)*S-DCON)/
     1        ((8.0D0/3.0D0)*C**4)
  550     CONTINUE
          PPF=DSQRT(DNU)*S/C
        ELSEIF(NU.EQ.6)THEN
C
C         AUGMENT THE RESULTS FOR THE NU = 6 CASE
C
          DCON=2.0D0*(DP-0.5D0)
          DARG=DPPF/DSQRT(DNU)
          Z=DATAN(DARG)
          DO650IPASS=1,MAXIT
            S=DSIN(Z)
            C=DCOS(Z)
            Z=Z-((1.0D0+0.5D0*C*C+0.375D0*C**4)*S-DCON)/
     1        ((15.0D0/8.0D0)*C**5)
  650     CONTINUE
          PPF=DSQRT(DNU)*S/C
        ENDIF
        GOTO9000
      ENDIF
C
C     CASE FOR FRACTIONAL DEGREES OF FREEDOM.  USE BISECTION
C     METHOD TO NUMERICALLY INVERT CDF FUNCTION.
C
 8000 CONTINUE
C
C     STEP 1: DETERMINE A BRACKETING INTERVAL.  USE 0 AS
C             EITHER THE LOWER OR UPPER LIMIT.
C
      MAXIT=3000
      IF(P.GT.0.5)THEN
        CALL NORPPF(P,XLTEMP,IPR)
CCCCC   CALL CAUPPF(P,XRTEMP)
        ARG=PI*P
        XRTEMP=-COS(ARG)/SIN(ARG)
      ELSEIF(P.LT.0.5)THEN
        CALL NORPPF(P,XRTEMP,IPR)
CCCCC   CALL CAUPPF(P,XLTEMP)
        ARG=PI*P
        XLTEMP=-COS(ARG)/SIN(ARG)
      ENDIF
      XL=DBLE(XLTEMP)
      XR=DBLE(XRTEMP)
      DNU=DBLE(ANU)
C
C  BISECTION METHOD
C
      DP=DBLE(P)
      IC = 0
      FXL = -DP
      FXR = 1.0D0 - DP
C
  105 CONTINUE
        X = (XL+XR)*0.5D0
        TERM1=1.0D0/(1.0D0 + X*X/DNU)
        TERM2=DNU/2.0D0
        TERM3=0.5D0
        TERM4=DBETAI(TERM1,TERM2,TERM3,IPR)
        IF(X.EQ.0.0D0)THEN
          DCDF=0.5D0
        ELSEIF(X.LE.0.0D0)THEN
          DCDF=0.5D0*TERM4
        ELSE
          DCDF=1.0D0 - 0.5D0*TERM4
        ENDIF
        P1=DCDF
        PPF=REAL(X)
        FCS = P1 - DP
        IF(FCS*FXL.GT.ZERO)THEN
          XL = X
          FXL = FCS
        ELSE
          XR = X
          FXR = FCS 
        ENDIF
        XRML = XR - XL
        IF(XRML.LE.SIG .AND. ABS(FCS).LE.EPS)GOTO9000
        IC = IC + 1
        IF(IC.LE.MAXIT)THEN
          GOTO105
        ELSE
          WRITE(IPR,130)
  130     FORMAT('***** ERROR--TPPF ROUTINE DID NOT CONVERGE')
        ENDIF
      GOTO9000
C
 9000 CONTINUE
      RETURN
      END
      SUBROUTINE VAR(X,N,IPR,XVAR,IDEBUG,IERROR)
C
C     PURPOSE--THIS SUBROUTINE COMPUTES THE
C              SAMPLE VARIANCE (WITH DENOMINATOR N-1)
C              OF THE DATA IN THE INPUT VECTOR X.
C              THE SAMPLE VARIANCE = (THE SUM OF THE
C              SQUARED DEVIATIONS ABOUT THE SAMPLE MEAN)/(N-1).
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 ARGUMENTS--XVAR   = THE SINGLE PRECISION VALUE OF THE
C                                COMPUTED SAMPLE VARIANCE.
C     OUTPUT--THE COMPUTED SINGLE PRECISION VALUE OF THE
C             SAMPLE VARIANCE (WITH DENOMINATOR N-1).
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--DOUBLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     REFERENCES--SNEDECOR AND COCHRAN, STATISTICAL METHODS,
C                 EDITION 6, 1967, PAGE 44.
C               --DIXON AND MASSEY, INTRODUCTION TO STATISTICAL
C                 ANALYSIS, EDITION 2, 1957, PAGE 38.
C               --MOOD AND GRABLE, 'INTRODUCTION TO THE THEORY
C                 OF STATISTICS, EDITION 2, 1963, PAGE 171.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 INFORMATION TECHNOLOGY LABORATORY
C                 NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY
C                 GAITHERSBURG, MD 20899-8980
C                 PHONE--301-975-2855
C     LANGUAGE--ANSI FORTRAN (1977)
C     VERSION NUMBER--82.6
C     ORIGINAL VERSION--JUNE      1972.
C     UPDATED         --SEPTEMBER 1975.
C     UPDATED         --NOVEMBER  1975.
C     UPDATED         --JUNE      1979.
C     UPDATED         --AUGUST    1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
      CHARACTER*4 IDEBUG
      CHARACTER*4 IERROR
C
      CHARACTER*4 ISUBN1
      CHARACTER*4 ISUBN2
C
C---------------------------------------------------------------------
C
      DOUBLE PRECISION DN
      DOUBLE PRECISION DX
      DOUBLE PRECISION DSUM
      DOUBLE PRECISION DMEAN
      DOUBLE PRECISION DVAR
C
      DIMENSION X(*)
C
C-----START POINT-----------------------------------------------------
C
      ISUBN1='VAR '
      ISUBN2='    '
      IERROR='NO'
      DMEAN=0.0
C
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
  999   FORMAT(1X)
        WRITE(IPR,51)
   51   FORMAT('***** AT THE BEGINNING OF VAR--')
        WRITE(IPR,52)IDEBUG,N
   52   FORMAT('IDEBUG,N = ',A4,2X,I8)
        DO55I=1,N
          WRITE(IPR,56)I,X(I)
   56     FORMAT('I,X(I) = ',I8,G15.7)
   55   CONTINUE
      ENDIF
C
C               ************************
C               **  COMPUTE VARIANCE  **
C               ************************
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      AN=N
C
      IF(N.LT.1)THEN
        IERROR='YES'
        WRITE(IPR,999)
        WRITE(IPR,111)
  111   FORMAT('***** ERROR IN VAR--')
        WRITE(IPR,112)
  112   FORMAT('      THE NUMBER OF OBSERVATIONS FOR THE RESPONSE ',
     1         'VARIABLE IS LESS THAN 1.')
        WRITE(IPR,117)N
  117   FORMAT('      THE NUMBER OF OBSERVATIONS = ',I8,'.')
        GOTO9000
      ELSEIF(N.EQ.1)THEN
        XVAR=0.0
        GOTO9000
      ENDIF
C
      HOLD=X(1)
      DO135I=2,N
        IF(X(I).NE.HOLD)GOTO139
  135 CONTINUE
      WRITE(IPR,999)
      WRITE(IPR,136)
  136 FORMAT('***** WARNING IN VAR--')
      WRITE(IPR,137)HOLD
  137 FORMAT('      THE RESPONSE VARIABLE HAS ALL ELEMENTS = ',G15.7)
      XVAR=0.0
      GOTO9000
  139 CONTINUE
C
C               *****************************
C               **  STEP 2--               **
C               **  COMPUTE THE VARIANCE.  **
C               *****************************
C
      DN=N
      DSUM=0.0D0
      DO200I=1,N
        DSUM=DSUM+DBLE(X(I))
  200 CONTINUE
      DMEAN=DSUM/DN
C
      DSUM=0.0D0
      DO300I=1,N
        DX=DBLE(X(I))
        DSUM=DSUM+(DX-DMEAN)**2
  300 CONTINUE
      DVAR=DSUM/(DN-1.0D0)
      XVAR=REAL(DVAR)
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT.      **
C               *****************
C
 9000 CONTINUE
      IF(IDEBUG.EQ.'ON')THEN
        WRITE(IPR,999)
        WRITE(IPR,9011)
 9011   FORMAT('***** AT THE END       OF VAR--')
        WRITE(IPR,9012)IERROR,DMEAN,XVAR
 9012   FORMAT('IERROR,DMEAN,XVAR = ',A4,2X,2G15.7)
      ENDIF
C
      RETURN
      END
