probfun.f [docs/17c/peakfqSA-jfe-distrib/src] Revision: default  Date:
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     SET OF SUBROUTINES TO DEAL WITH PROBABILITY FUNCTIONS
C
C     COPYRIGHT, TIMOTHY A. COHN, 1995
C     PROPERTY OF US GOVERNMENT, GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C     
C
C     AUTHOR.......TIM COHN
C           DATE.....JANUARY 5, 1995
C
C     Minor modifications and additions
C     by John England, USACE john.f.england@usace.army.mil
C
C     modifications by John England denoted with C JFE
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C    07-SEP-2017 extended H-S routine arrays to 20000
C          
C
C
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
CC
C     N.B.  IMSL SUBROUTINES ARE USED:
C                 DGAMMA
C                 DGAMDF
C                 DNORDF
C                 DNORIN
C                 DCHIIN
C
C     N.B.  CDFLIB ROUTINES ARE USED
C                 CDFBET
C
C     ROUTINES AVAILABLE
C           STANDARD NORMAL DISTRIBUTION
C                 FP_Z_PDF(X)
C                 FP_Z_PDF_LN(X)
C                 FP_Z_PDF_TRA(X,T)
C                 FP_Z_CDF(X)
C                 FP_Z_CDF_TRA(X,T)
C                 FP_Z_ICDF(P)
C                 FP_Z_ICDF_TRA(P,T)
C                 FP_Z_MOM(K)
C                 FP_Z_MOM_TRA(K,T)
C                 FP_Z_MOM_TRB(K,T)
C                 SP_Z(X,F,PHI,LOGPHI,A,G)
C
C           2 PARAMETER NORMAL
C                 FP_N2_PDF(X,PARMS)
C                 FP_N2_CDF(X,PARMS)
C                 FP_N2_ICDF(P,PARMS)
C
C           GAMMA DISTRIBUTION
C                 FP_G1_ICDF(P,ALPHA)
C                 FP_G2_ICDF(P,PARMS)  -- PARMS=(ALPHA,BETA)
C                 FP_G3_ICDF(P,PARMS)  -- PARMS=(TAU,ALPHA,BETA)
C                 FP_G1_MOM_TRA(ALPHA,T,K) 
C                 FP_G2_MOM_TRA(PARMS,T,K) 
C                 FP_G3_MOM_TRA(PARMS,T,K) 
C
C           NON-CENTRAL STUDENTS T DISTRIBUTION
C                 FP_NCT_CDF(X,NU,DELTA)
C                 FP_NCT_ICDF(P,NU,DELTA)
C                 FP_T_CDF(T,DF,NCP)  ** PROBABLY OBSOLETE
C                 FP_T_ICDF(P,DF,NCP)  ** PROBABLY OBSOLETE
C
C           BETA DISTRIBUTION
C                 FP_B_CDF(X,ALPHA,BETA)
C                 FP_B_ICDF(P,ALPHA,BETA)
C
C     NAMING CONVENTIONS:
C
C     1)    TYPE OF ROUTINE
C                 FP_         =     DOUBLE PRECISION FUNCTIONS
C                 SP_         =     SUBROUTINES
C
C     2)    NEXT COMES LETTERS DEFINING THE DISTRIBUTION:
C
C                 Z_          =     PHI(X)                              STANDARD NORMAL DISTRIBUTION
C                 N_          =     N(MU,SIGMA)                   2-P NORMAL DISTRIBUTION 
C                 G1_         =     GAMMA(X,ALPHA)                1-P GAMMA DISTRIBUTION 
C                 G2_         =     GAMMA(X,ALPHA,BETA)           2-P GAMMA DISTRIBUTION 
C                 G3_         =     GAMMA(X,ALPHA,BETA,TAU) 3-P GAMMA DISTRIBUTION 
C
C     3)    NEXT COMES FUNCTION TYPE:
C                 
C                 PDF_        =     PROBABILITY DENSITY FUNCTION
C                 CDF_        =     CUMULATIVE DENSITY FUNCTION
C                 ICDF_       =     INVERSE CDF
C                 MOM_        =     MOMENTS
C                 HAZ_        =     F/(1-F) 
C                 RAT_        =     F/F
C
C     4)    NEXT COMES WHETHER LOGS ARE USED
C
C                 LN_               =     USE NATURAL LOG
C
C     5)    NEXT COMES FUNCTION TYPE
C
C                 JAC_  =     JACOBIAN (VECTOR OF FIRST DERIVATIVES WRT PARAMS)
C                 HESS_ =     HESSIAN OF PDF (MATRIX OF SECOND DERIVATIVES WRT PARAMS)
C                 
C     6)    NEXT COMES TRUNCATION
C
C                 TRA_  =     DISTRIBUTION TRUNCATED ABOVE (I.E. ALL VALUES "LESS THAN")
C                 TRB_  =     DISTRIBUTION TRUNCATED BELOW
C 
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_Z_PDF(X)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE PDF OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      CALL SP_Z(X,F,PHI,XLNPHI,RAT,HAZ)
      
      FP_Z_PDF    =     F
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_Z_CDF(X)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      CALL SP_Z(X,F,PHI,XLNPHI,RAT,HAZ)
      
      FP_Z_CDF    =     PHI
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_Z_ICDF(P)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE CDF OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      FP_Z_ICDF   =     DNORIN(P)
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_Z_MOM(K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE MOMENTS OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION RES(1000)
      DATA RES/0,1,0,3,0,15,0,105,0,991*-99/
      
      IF(K .GT. 9) THEN
        DO 10 I=1,K
          RES(K) = (I-1.D0)*RES(K-2)
10    CONTINUE
      ENDIF

      FP_Z_MOM    =     RES(K)
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_N2_PDF(X,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE PDF OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION PARMS(2),MU,SIGMA

            MU          =     PARMS(1)
            SIGMA =     PARMS(2)    
            Z     =     (X-MU)/SIGMA
      CALL SP_Z(Z,F,PHI,XLNPHI,RAT,HAZ)
      
      FP_N2_PDF   =     F/SIGMA
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_N2_CDF(X,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION PARMS(2),MU,SIGMA
      
            MU          =     PARMS(1)
            SIGMA =     PARMS(2)
            Z     =     (X-MU)/SIGMA
      
      CALL SP_Z(Z,F,PHI,XLNPHI,RAT,HAZ)
      
      FP_N2_CDF   =     PHI
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_N2_ICDF(P,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE CDF OF STANDARD NORMAL DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION PARMS(2),MU,SIGMA

            MU          =     PARMS(1)
            SIGMA =     PARMS(2)
      
      FP_N2_ICDF  =     MU + SIGMA * DNORIN(P)
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G1_PDF(X,ALPHA)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE PDF OF 1-P GAMMA DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      IF(X .GT. 0.D0 .AND. ALPHA .GT. 0.D0) THEN
            IF((ALPHA+X) .LT. 100.D0) THEN
                  FP_G1_PDF   =     X**(ALPHA-1.D0) * 
     1                                 EXP(-X)/DGAMMA(ALPHA)
            ELSE
                  TEMP        =     (ALPHA-1.D0) * LOG(X) - X -  
     1                                 DLNGAM(ALPHA)
                  FP_G1_PDF   =     EXP(TEMP)
            ENDIF 
      ELSE
            FP_G1_PDF   =     0.D0
      ENDIF
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G2_PDF(X,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE PDF OF 2-P GAMMA DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION PARMS(2)
      
      ALPHA =     PARMS(1)
      BETA  =     PARMS(2)
      
      FP_G2_PDF   =     FP_G1_PDF(X/BETA,ALPHA)/ABS(BETA)
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G3_PDF(X,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE PDF OF 3-P GAMMA DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION PARMS(3)
      
      TAU         =     PARMS(1)
      ALPHA =     PARMS(2)
      BETA  =     PARMS(3)
      
      FP_G3_PDF   =     FP_G2_PDF(X-TAU,PARMS(2))
      
      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~

      DOUBLE PRECISION FUNCTION FP_G1_CDF(X,ALPHA)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF 1-P GAMMA DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      FP_G1_CDF   =     DGAMDF(X,ALPHA)
      
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~

      DOUBLE PRECISION FUNCTION FP_G2_CDF(X,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF 2-P GAMMA DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION PARMS(2)
      
      ALPHA =     PARMS(1)
      BETA  =     PARMS(2)
      
      IF(BETA .GT. 0.D0) THEN
            FP_G2_CDF   =     FP_G1_CDF(X/BETA,ALPHA)
      ELSE IF(BETA .LT. 0.D0) THEN
            FP_G2_CDF   =     1.D0 - FP_G1_CDF(X/BETA,ALPHA)
      ELSE
        IF(X .GT. 0.D0) THEN
              FP_G2_CDF   =     1.D0
        ELSE
              FP_G2_CDF   =     0.D0
      ENDIF
      ENDIF
      
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~

      DOUBLE PRECISION FUNCTION FP_G3_CDF(X,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF 2-P GAMMA DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION PARMS(3)
      
      TAU         =     PARMS(1)
C     ALPHA =     PARMS(2)
C     BETA  =     PARMS(3)
      
            FP_G3_CDF   =     FP_G2_CDF(X-TAU,PARMS(2))

      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G1_ICDF(P,ALPHA)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE OF GAMMA(ALPHA)
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      IF(P .GT. 0.999999999) THEN
            FP_G1_ICDF  =     9.D99
            RETURN
      ELSE IF(P .LT. 0.000000001) THEN
            FP_G1_ICDF  =     0.D0
            RETURN      
      ENDIF
      
            FP_G1_ICDF  =     0.5*DCHIIN(P,2.D0*ALPHA)
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C
      DOUBLE PRECISION FUNCTION FP_G2_ICDF(P,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE OF GAMMA(P,ALPHA,BETA)
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION PARMS(2)
      
      ALPHA =     PARMS(1)
      BETA  =     PARMS(2)
      
      IF(BETA .GT. 0.D0) THEN
            FP_G2_ICDF  =     BETA * FP_G1_ICDF(P,ALPHA)
      ELSE IF(BETA .LT. 0.D0) THEN
            FP_G2_ICDF  =     BETA * FP_G1_ICDF(1.D0-P,ALPHA)
      ELSE
            WRITE(*,*) ' INVALID PARAMETERS;  BETA = 0'
            RETURN
      ENDIF

      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G3_ICDF(P,PARMS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE OF GAMMA(P,TAU,ALPHA,BETA)
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DIMENSION PARMS(3)
      
      TAU         =     PARMS(1)
C     ALPHA =     PARMS(2)
C     BETA  =     PARMS(3)
      
      FP_G3_ICDF  =     TAU + FP_G2_ICDF(P,PARMS(2))
      
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G1_MOM_TRA(ALPHA,T,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE THE EXPECTED VALUE OF THE K-TH MOMENT OF 
C     A 1-PARAMETER GAMMA VARIATE CENSORED ABOVE AT T
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C

      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE

       T1   =     MAX(0.D0,T)
       ANS  =     DGAMDF(T1,ALPHA+K)/MAX(1.D-50,DGAMDF(T1,ALPHA))

C     THIS SECTION COMPUTES RATIO GAMMA(ALPHA+K)/GAMMA(ALPHA)

      DO 10 J=0,K-1
            ANS   =     ANS*(ALPHA+J)
10    CONTINUE

      FP_G1_MOM_TRA     =     ANS

      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G2_MOM_TRA(PARMS,T,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE THE EXPECTED VALUE OF THE K-TH MOMENT OF 
C     A 2-PARAMETER GAMMA VARIATE CENSORED ABOVE AT T
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C

      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(2)
      
      ALPHA =     PARMS(1)
      BETA  =     PARMS(2)
      
      IF(BETA .EQ. 0.D0) THEN
            WRITE(*,*) ' INVALID PARAMETERS;  BETA = 0'
            RETURN
      ELSE IF(BETA .GT. 0.D0) THEN
            FP_G2_MOM_TRA = BETA**K * FP_G1_MOM_TRA(ALPHA,T/BETA,K)
      ELSE
            FP_G2_MOM_TRA = BETA**K * FP_G1_MOM_TRB(ALPHA,T/BETA,K)
      ENDIF
            
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G3_MOM_TRA(PARMS,T,K)
C****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE KTH NON-CENTRAL MOMENT OF PEARSON TYPE III 
C     RANDOM VARIABLE CENSORED AT T WITH PARAMETERS PARMS
C
C
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(3)

            TAU         =     PARMS(1)
            ALPHA =     PARMS(2)
            BETA  =     PARMS(3)

            SUM   =     FP_G2_MOM_TRA(PARMS(2),T-TAU,K)
            
      IF(TAU .NE. 0.D0) THEN
            DO 10 J=0,K-1
                  SUM   =     SUM + CHOOSE(K,J)*TAU**(K-J)*
     1                       FP_G2_MOM_TRA(PARMS(2),T-TAU,J)
10          CONTINUE
      ENDIF
      
            FP_G3_MOM_TRA     =     SUM

      RETURN
      END

      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G1_MOM_TRB(ALPHA,T,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE THE EXPECTED VALUE OF THE K-TH MOMENT OF 
C     A 1-PARAMETER GAMMA VARIATE CENSORED BELOW AT T
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C

      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE

            T1    =     MAX(0.D0,T)
            ANS         =     (1.D0 - DGAMDF(T1,ALPHA+K))/
     1                       MAX(1.D-50, 1.D0 - DGAMDF(T1,ALPHA))

C     THIS SECTION COMPUTES RATIO GAMMA(ALPHA+K)/GAMMA(ALPHA)

      DO 10 J=0,K-1
            ANS   =     ANS*(ALPHA+J)
10    CONTINUE

      FP_G1_MOM_TRB     =     ANS

      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G2_MOM_TRB(PARMS,T,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE THE EXPECTED VALUE OF THE K-TH MOMENT OF 
C     A 2-PARAMETER GAMMA VARIATE CENSORED ABOVE AT T
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C

      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(2)
      
      ALPHA =     PARMS(1)
      BETA  =     PARMS(2)
      
      IF(BETA .EQ. 0.D0) THEN
            WRITE(*,*) ' INVALID PARAMETERS;  BETA = 0'
            RETURN
      ELSE IF(BETA .GT. 0.D0) THEN
            FP_G2_MOM_TRB = BETA**K * FP_G1_MOM_TRB(ALPHA,T/BETA,K)
      ELSE
            FP_G2_MOM_TRB = BETA**K * FP_G1_MOM_TRA(ALPHA,T/BETA,K)
      ENDIF
            
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G3_MOM_TRB(PARMS,T,K)
C****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE KTH NON-CENTRAL MOMENT OF PEARSON TYPE III 
C     RANDOM VARIABLE CENSORED AT T WITH PARAMETERS PARMS
C
C
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(3)

            TAU         =     PARMS(1)
            ALPHA =     PARMS(2)
            BETA  =     PARMS(3)

            SUM   =     FP_G2_MOM_TRB(PARMS(2),T-TAU,K)
            
      IF(TAU .NE. 0.D0) THEN
            DO 10 J=0,K-1
                  SUM   =     SUM + CHOOSE(K,J)*TAU**(K-J)*
     1                       FP_G2_MOM_TRB(PARMS(2),T-TAU,J)
10          CONTINUE
      ENDIF
      
            FP_G3_MOM_TRB     =     SUM

      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G1_MOM_TRC(ALPHA,TL,TU,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE THE EXPECTED VALUE OF THE K-TH MOMENT OF 
C     A 1-PARAMETER GAMMA VARIATE BETWEEN TL AND TU WITH PARAMETER ALPHA
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C     
C     MODIFIED...DECEMBER 23, 1996 (TAC)
C     MODIFIED...JULY 25, 1998 (TAC)
C

      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE

            TL1   =     MAX(MIN(0.D0,TU),TL)
      IF(TL1 .EQ. TU) THEN
          ANS = TL1**K
      ELSE
        DOWN    =  DGAMDF(TU,ALPHA) - DGAMDF(TL1,ALPHA)
      IF(DOWN .GT. 0.D0) THEN
        UP    =  DGAMDF(TU,ALPHA+K) - DGAMDF(TL1,ALPHA+K)
        ANS    =  UP/DOWN
C     THIS SECTION COMPUTES RATIO GAMMA(ALPHA+K)/GAMMA(ALPHA)
            DO 10 J=0,K-1
                  ANS   =     ANS*(ALPHA+J)
10          CONTINUE
      ELSE
        IF(TL1 .GT. ALPHA) THEN
           ANS = TL1**K
        ELSE
           ANS = TU**K
        ENDIF
      ENDIF
      ENDIF

      FP_G1_MOM_TRC     =     ANS

      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G2_MOM_TRC(PARMS,TL,TU,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE THE EXPECTED VALUE OF THE K-TH MOMENT OF 
C     A 2-PARAMETER GAMMA VARIATE BETWEEN TL AND TU WITH PARAMETERS PARMS

C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 28, 1994
C

      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(2)
      
      ALPHA =     PARMS(1)
      BETA  =     PARMS(2)
      
      IF(BETA .EQ. 0.D0) THEN
            WRITE(*,*) ' INVALID PARAMETERS;  BETA = 0'
            RETURN
      ELSE IF(BETA .GT. 0.D0) THEN
            FP_G2_MOM_TRC     =     BETA**K *  
     1            FP_G1_MOM_TRC(ALPHA,TL/BETA,TU/BETA,K)
      ELSE
            FP_G2_MOM_TRC     =     BETA**K *  
     1            FP_G1_MOM_TRC(ALPHA,TU/BETA,TL/BETA,K)
      ENDIF
            
      RETURN
      END
      
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G3_MOM_TRC(PARMS,TL,TU,K)
C****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE KTH NON-CENTRAL MOMENT OF PEARSON TYPE III 
C     RANDOM VARIABLE BETWEEN TL AND TU WITH PARAMETERS PARMS
C
C
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(3)

            TAU         =     PARMS(1)
            ALPHA =     PARMS(2)
            BETA  =     PARMS(3)

            SUM   =     FP_G2_MOM_TRC(PARMS(2),TL-TAU,TU-TAU,K)
            
      IF(TAU .NE. 0.D0) THEN
            DO 10 J=0,K-1
                  SUM   =     SUM + CHOOSE(K,J)*TAU**(K-J)*
     1                       FP_G2_MOM_TRC(PARMS(2),TL-TAU,TU-TAU,J)
10          CONTINUE
      ENDIF
      
            FP_G3_MOM_TRC     =     SUM

      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_G3_MOM_TRCO(PARMS,X_OFF,TL,TU,K)
C****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE KTH MOMENT, CENTRALIZED AROUND X_OFF, OF PEARSON TYPE III 
C     RANDOM VARIABLE BETWEEN TL AND TU WITH PARAMETERS PARMS
C
C
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
C      SAVE
      DIMENSION PARMS(3)

            TAU   =     PARMS(1)
            ALPHA =     PARMS(2)
            BETA  =     PARMS(3)

            SUM   =     FP_G2_MOM_TRC(PARMS(2),TL-TAU,TU-TAU,K)
            
      IF(TAU - X_OFF .NE. 0.D0) THEN
            DO 10 J=0,K-1
                  SUM   =     SUM + CHOOSE(K,J)*(TAU-X_OFF)**(K-J)*
     1                       FP_G2_MOM_TRC(PARMS(2),TL-TAU,TU-TAU,J)
10          CONTINUE
      ENDIF
      
            FP_G3_MOM_TRCO    =     SUM

      RETURN
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
       SUBROUTINE SP_Z(XSI,F,PHI,LOGPHI,A,G)
C======================================================================
C
C      SUBROUTINE TO COMPUTE GAUSSIAN FUNCTIONS AND THEIR
C      DERIVATIVES, WHERE XSI IS A STANDARD NORMAL VARIATE
C
C           PHI       =  STANDARD NORMAL CDF
C           LOGPHIT   =  LOG(PHIT)
C           AT        =  F/PHIT
C           GT        =  F/(1.D0-PHIT)
C
C           THE ROUTINE WORKS ACCURATELY OVER A VERY LONG RANGE, AND PROVIDES
C           A 2ND ORDER CONTINUITY CORRECTION WHERE ESTIMATION METHODS CHANGE.
C
C
C        AUTHOR......TIM COHN
C        DATE........APRIL 1, 1987
C          ERROR CORRECTED OCTOBER 4, 1992 (TAC)
C                 IMPROVED VERSION....DECEMBER 20, 1994
C
C======================================================================
       IMPLICIT DOUBLE PRECISION (A-Z)
C      SAVE

       DATA LIMIT/6.5D0/,LIMITI/6.D0/,CONST/0.3989422804014327/
         
           Z(B)      =  1.D0/( SQRT(B)*(1.D0-1.D0*B*(1.D0-3.D0*B*
     1                  (1.D0-5.D0*B*(1.D0-7.D0*B*(1.D0-9.D0*B*
     2                  (1.D0-11.D0*B*(1.D0-13.D0*B*(1.D0-15.D0*B*
     3                  (1.D0-17.D0*B*(1.D0-19.D0*B*(1.D0-21.D0*B*
     4                  (1.D0-23.D0*B*(1.D0-12.5D0*B))))))))))))) )
              
         F         =  CONST*EXP(-0.5*XSI**2)
       IF(XSI .LT. -LIMITI) THEN
         PHI       =  F/Z(1.D0/XSI**2) 
         A         =  Z(1.D0/XSI**2)
         LOGPHI    =  LOG(CONST/A) - 0.5*XSI**2
         G         =  F/(1.D0-PHI) 
       ELSE IF (XSI .GT. LIMITI) THEN
         PHI       =  1.D0 - F/Z(1.D0/XSI**2) 
         LOGPHI    =  LOG(PHI) 
         A         =  F/PHI 
         G         =  Z(1.D0/XSI**2)
       ENDIF
         
       IF(ABS(XSI) .LT. LIMIT) THEN
           PHIT       =  DNORDF(XSI)
           LOGPHIT    =  LOG(PHIT)
           AT         =  F/PHIT
           GT         =  F/(1.D0-PHIT)
             IF(ABS(XSI) .GT. LIMITI) THEN
                  IF(ABS(XSI) .GT. LIMIT) THEN
                        ARG   =     1.D0
                  ELSE
                     ARGT      =  (ABS(XSI)-(LIMIT+LIMITI)/2.D0) 
     1                                   /(LIMIT-LIMITI)
                     ARG       =  (1.D0+SIN(3.14159*ARGT))/2.D0
                  ENDIF
             ELSE
               ARG       =  0.D0
             ENDIF
               PHI       =  ARG*PHI + (1.D0-ARG)*PHIT
               LOGPHI    =  ARG*LOGPHI + (1.D0-ARG)*LOGPHIT
               A         =  ARG*A + (1.D0-ARG)*AT
               G         =  ARG*G + (1.D0-ARG)*GT
       ENDIF
       RETURN
       END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION CHOOSE(N,K)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE COMBINATORIAL CHOOSE FUNCTION
C
C                 / N \
C                 \ K /
C
C     N.B.  N AND K MUST BE INTEGERS;  CHOOSE IS DOUBLE PRECISION
C
C     AUTHOR.....TIM COHN
C     DATE.......DECEMBER 21, 1994
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      INTEGER N,K,TAB1(0:12,0:12)
      
      DATA TAB1/
     * 1,  12*0,
     * 1,1,  11*0,
     * 1,2,1,  10*0,
     * 1,3,3,1,  9*0,
     * 1,4,6,4,1,  8*0,
     * 1,5,10,10,5,1,  7*0,
     * 1,6,15,20,15,6,1,  6*0,
     * 1,7,21,35,35,21,7,1,  5*0,
     * 1,8,28,56,70,56,28,8,1,  4*0,
     * 1,9,36,84,126,126,84,36,9,1,  3*0,
     * 1,10,45,120,210,252,210,120,45,10,1,  2*0,
     * 1,11,55,165,330,462,462,330,165,55,11,1,  1*0,
     * 1,12,66,220,495,792,924,792,495,220,66,12,1 /
      
      IF(N .LT. MAX(0,K) .OR. K .LT. 0) THEN
            WRITE(*,*) ' ERROR IN CHOOSE: (N,K) ',N,K
            STOP
      ENDIF
      
      IF(N .LE. 12) THEN
            CHOOSE      =     TAB1(K,N)
      ELSE IF(N .LT. 170.D0) THEN
            CHOOSE      =     DFAC(N)/(DFAC(N-K)*DFAC(K))
      ELSE
            ARG         =     DLNGAM(N+1.D0)-DLNGAM(N-K+1.D0)- 
     1                                   DLNGAM(K+1.D0)
            IF(ARG .LT. 500.D0) THEN
                  CHOOSE      =     EXP(ARG)
            ELSE
                  CHOOSE      =     -99.D0
            ENDIF
      ENDIF
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      SUBROUTINE SP_G1_CDF_HESS(X,ALPHA_IN,GRAD,HESS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE HESSIAN AND GRADIENT OF 1-P GAMMA D''N
C
C     AUTHOR.....TIM COHN
C     DATE.......JANUARY 30, 1995
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     ORDER OF PARAMETERS IS (X, ALPHA)
C
C           GRAD(1)     =     D_F/D_X
C           GRAD(2)     =     D_F/D_ALPHA
C
C           HESS(1,1)   =     D2_F/D_X2
C           HESS(1,2)   =     HESS(2,1)   =     D2_F/D_X D_ALPHA 
C           HESS(2,2)   =     D2_F/D_ALPHA2
C
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      EXTERNAL FCNGAM
      
      DIMENSION GRAD(2),HESS(2,2)
      
      COMMON /TACGAM001/ALPHA,M
      
      DATA ERRABS/0.D0/,ERREL/1.D-5/
      
      ALPHA =     ALPHA_IN
      
            P     =     FP_G1_CDF(X,ALPHA)
            
            FPDF        =     FP_G1_PDF(X,ALPHA)
            GRAD(1)           =     FPDF
            HESS(1,1)   =     FPDF * ( (ALPHA-1.D0)/X - 1.D0)
            HESS(1,2)   =     FPDF * ( LOG(X) - DPSI(ALPHA) )
            HESS(2,1)   =     HESS(1,2)
C
      M     =     1
            CALL DQDAGS(FCNGAM,0.D0,P,   ERRABS,ERREL,FP,ERREST)
            CALL DQDAGS(FCNGAM,0.D0,1.D0,ERRABS,ERREL,GP,ERREST)
      M     =     2
            CALL DQDAGS(FCNGAM,0.D0,P,   ERRABS,ERREL,FPP,ERREST)
            CALL DQDAGS(FCNGAM,0.D0,1.D0,ERRABS,ERREL,GPP,ERREST)
C
C
C     N.B.  P IS F
C           BECAUSE G   =     1.D0
C                 D1    =     FP*G - GP*P
C                 D2    =     ( ( (FPP * G + GP * FP) ) * G**2 - 2 * G * GP * (FP*G - GP*P) )/G**4
C
      GRAD(2)     =     FP - GP*P

      HESS(2,2)   = ( ( (FPP +     GP * FP) - (GPP * P + GP * FP) )    
     1                  - 2 *     GP * (FP   - GP*P) )
      
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION FCNGAM(P)
C      SAVE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      
      COMMON /TACGAM001/ALPHA,M
      
      X           =     FP_G1_ICDF(P,ALPHA)
      FCNGAM      =     LOG(X)**M
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION SIGN1(ARG)
C****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C
      DOUBLE PRECISION ARG
C      SAVE

      IF(ARG .LT. 0.D0) THEN
            SIGN1 =     -1.D0
      ELSE
            SIGN1 =      1.D0
      ENDIF
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      SUBROUTINE SP_G1_CDF_LN_HESS(X,ALPHA,GRAD,HESS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE HESSIAN AND GRADIENT OF LOGARITHM OF 1-P GAMMA D''N
C
C     AUTHOR.....TIM COHN
C     DATE.......MARCH 1, 1995
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     ORDER OF PARAMETERS IS (X, ALPHA)
C
C           GRAD(1)     =     D_LOG(F)/D_X
C           GRAD(2)     =     D_LOG(F)/D_ALPHA
C
C           HESS(1,1)   =     D2_LOG(F)/D_X2
C           HESS(1,2)   =     HESS(2,1)   =     D2_LOG(F)/D_X D_ALPHA 
C           HESS(2,2)   =     D2_LOG(F)/D_ALPHA2
C
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      DIMENSION GRAD(2),HESS(2,2),G1(2),H2(2,2)
      
      F     =     FP_G1_CDF(X,ALPHA)
      
      IF(F .GT. 0.D0) THEN
                  CALL SP_G1_CDF_HESS(X,ALPHA,G1,H2)
            DO 10 I=1,2
                  GRAD(I)           =     G1(I)/F
            DO 10 J=1,2
                  HESS(I,J)   =     ( H2(I,J)*F - G1(I)*G1(J) )/F**2
10          CONTINUE
      ELSE
            DO 20 I=1,2
                  GRAD(I)           =     0.D0
            DO 20 J=1,2
                  HESS(I,J)   =     0.D0
20          CONTINUE
      ENDIF
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      SUBROUTINE SP_G3_CDF_LN_HESS(X,PARMS,GRAD,HESS)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE HESSIAN AND GRADIENT OF LOGARITHM OF 3-P GAMMA D''N
C
C     AUTHOR.....TIM COHN
C     DATE.......MARCH 1, 1995
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     ORDER OF PARAMETERS IS (TAU, ALPHA, BETA)
C     HESS IS A SYMMETRIC MATRIX
C
C           GRAD(1)     =     D_LOG(F)/D_TAU
C           GRAD(2)     =     D_LOG(F)/D_ALPHA
C           GRAD(3)     =     D_LOG(F)/D_BETA
C
C           HESS(1,1)   =     D2_LOG(F)/D_TAU^2
C           HESS(1,2)   =     D2_LOG(F)/D_TAU D_ALPHA 
C           HESS(1,3)   =     D2_LOG(F)/D_TAU D_BETA 
C           HESS(2,2)   =     D2_LOG(F)/D_ALPHA^2
C           HESS(2,3)   =     D2_LOG(F)/D_ALPHA D_BETA 
C           HESS(3,3)   =     D2_LOG(F)/D_BETA^2
C
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      DIMENSION PARMS(3),GRAD(3),HESS(3,3),G1(2),H2(2,2)
      
            TAU         =     PARMS(1)
            ALPHA =     PARMS(2)
            BETA  =     PARMS(3)
            
            Y           =     (X-TAU)/BETA
                        
      CALL SP_G1_CDF_LN_HESS(Y,ALPHA,G1,H2)
      
                  GRAD(1)           =     G1(1) * (-1.D0/BETA)
                  GRAD(2)           =     G1(2)
                  GRAD(3)           =     G1(1) * (-Y/BETA)
                  
                  HESS(1,1)   =     H2(1,1) * (1.D0/BETA**2)
                  HESS(1,2)   =     H2(1,2) * (-1.D0/BETA)
                  HESS(1,3)   =     H2(1,1) * (Y/BETA**2) + G1(1) *  
     1                                   (1.D0/BETA**2)
                  HESS(2,2)   =     H2(2,2)
                  HESS(2,3)   =     H2(1,2) * (-Y/BETA)
                  HESS(3,3)   =     H2(1,1) * (-Y/BETA)**2 + G1(1) *  
     1                                   (2*Y/BETA**2)
                  
      DO 10 I=1,3
            DO 10 J=I+1,3
                  HESS(J,I)   =     HESS(I,J)
10    CONTINUE

      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      SUBROUTINE DPDM(PARMS,JACT)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE JACOBIAN DP/DM TRANSFORMATION CONVERTING
C        PEARSON TYPE 3 NON-CENTRAL MOMENTS (M1,M2,M3)
C        INTO PARAMETERS (TAU,ALPHA,BETA) 
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 10, 1999
C
C    N.B.  JAC IS COMPUTED;  JACT = TRANSPOSE(JAC) IS RETURNED
C
C            JAC(1,1) =  DT/DM1
C            JAC(1,2) =  DA/DM1
C            JAC(1,3) =  DB/DM1
C
C            JAC(2,1) =  DT/DM2
C            JAC(2,2) =  DA/DM2
C            JAC(2,3) =  DB/DM2
C
C            JAC(3,1) =  DT/DM3
C            JAC(3,2) =  DA/DM3
C            JAC(3,3) =  DB/DM3
C
C    N.B.: THIS CODE APPEARS TO HANDLE ALL CASES CORRECTLY (TAC 4/10/99)
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      
      DOUBLE PRECISION PARMS(3),M(3),M1,M2,M3,S,A,B,D,T
      DOUBLE PRECISION JAC(3,3),JACT(3,3)
            
      CALL P2MN(PARMS,M)
        M1  =  M(1)
        M2  =  M(2)
        M3  =  M(3)
      
C   V IS THE VARIANCE

      V        =  M2 - M1**2
        DVDM1  =  -2.D0*M1
        DVDM2  =  1.D0
        DVDM3  =  0.D0

C   S  IS THE SKEW       
      S        =  M3 - 3.D0*M2*M1 + 2.D0*M1**3
        DSDM1  =  -3.D0*M2 + 6.D0*M1**2
        DSDM2  =  -3.D0*M1
        DSDM3  =  1.D0

C   A  IS THE PARAMETER

      A        =  4.D0*(V**3)/S**2
        DADM1  =  4.D0*(3.D0*V**2*DVDM1 * S**2 - 2*S*DSDM1*V**3)/S**4
        DADM2  =  4.D0*(3.D0*V**2*DVDM2 * S**2 - 2*S*DSDM2*V**3)/S**4
        DADM3  =  4.D0*(3.D0*V**2*DVDM3 * S**2 - 2*S*DSDM3*V**3)/S**4
        
        
C   D IS B**2

      D        =  V/A
        DDDM1  =  (DVDM1*A-DADM1*V)/A**2
        DDDM2  =  (DVDM2*A-DADM2*V)/A**2
        DDDM3  =  (DVDM3*A-DADM3*V)/A**2
        
C   B IS SIGN(S)*SQRT(D)

      B        =  SIGN(1.D0,S)*SQRT(D)
        DBDM1  =  0.5*DDDM1/B
        DBDM2  =  0.5*DDDM2/B
        DBDM3  =  0.5*DDDM3/B

C   T

      T        =  M1 - A*B
        DTDM1  =  1.D0 - (DADM1*B + A*DBDM1)
        DTDM2  =       - (DADM2*B + A*DBDM2)
        DTDM3  =       - (DADM3*B + A*DBDM3)

      JAC(1,1) =  DTDM1
      JAC(1,2) =  DADM1
      JAC(1,3) =  DBDM1

      JAC(2,1) =  DTDM2
      JAC(2,2) =  DADM2
      JAC(2,3) =  DBDM2

      JAC(3,1) =  DTDM3
      JAC(3,2) =  DADM3
      JAC(3,3) =  DBDM3
      
      DO 10 I=1,3
        DO 10 J=1,3
          JACT(I,J) = JAC(J,I)
10    CONTINUE

      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C        
        SUBROUTINE P2M(P,M)
C        SAVE
        IMPLICIT NONE
        DOUBLE PRECISION P(3),M(3),T,A,B
        T=P(1)
        A=P(2)
        B=P(3)
        M(1) =  A*B + T
        M(2) =  A*B**2
        M(3) =  2.D0/SQRT(A)
        IF(B .LT. 0.D0) M(3) = -M(3)
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
        SUBROUTINE M2P(M,P)
        IMPLICIT NONE
C        SAVE
        DOUBLE PRECISION P(3),A,B,T,M(3),M1,M2,M3
        M1=M(1)
        M2=M(2)
        M3=M(3)
        A = 4.D0/M3**2
        B = SQRT(M2/A)
        IF(M3 .LT. 0.D0) B=-B
        T = M1 - A*B
        P(1)  =  T
        P(2)  =  A
        P(3)  =  B
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C        
        SUBROUTINE P2MN(P,M)
C        SAVE
        IMPLICIT NONE
        DOUBLE PRECISION P(3),M(3),T,A,B
        T=P(1)
        A=P(2)
        B=P(3)
        M(1) =  A*B + T
        M(2) =  A*(1 + A)*B**2 + 2*A*B*T + T**2
        M(3) =  A*(1 + A)*(2 + A)*B**3 + 3*A*(1 + A)*B**2*T + 
     1             3*A*B*T**2 + T**3
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
        SUBROUTINE MN2P(M,P)
        IMPLICIT NONE
C      SAVE
        DOUBLE PRECISION P(3),M(3),M1,M2,M3
        M1=M(1)
        M2=M(2)
        M3=M(3)
        P(1) = (M1**2*M2 - 2*M2**2 + M1*M3)/(2*M1**3 - 3*M1*M2 + M3)
        P(2) = -4*(M1**2 - M2)**3/(2*M1**3 - 3*M1*M2 + M3)**2
        P(3) = (2*M1**3 - 3*M1*M2 + M3)/(-2*M1**2 + 2*M2)
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C        
        SUBROUTINE M2MN(M,MN)
C      SAVE
        IMPLICIT NONE
        DOUBLE PRECISION M(3),MN(3)
        MN(1) =  M(1)
        MN(2) =  M(2)+M(1)**2
        MN(3) =  M(3)*M(2)**1.5 + 3.D0*MN(2)*M(1) - 2.D0 * M(1)**3
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C        
        SUBROUTINE MN2M(MN,M)
C      SAVE
        IMPLICIT NONE
        DOUBLE PRECISION M(3),MN(3)
        M(1) =  MN(1)
        M(2) =  MN(2)-MN(1)**2
        M(3) =  (MN(3) - 3.D0*MN(2)*MN(1) + 2.D0 * MN(1)**3)
     1              /M(2)**1.5
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION D_G1_INV(PQ,ALPHA)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C     PROGRAM TO COMPUTE DERIVATIVE WRT ALPHA OF 1-PARAMETER GAMMA QUANTILE
C         FUNCTION
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 10, 1999
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C           PQ        =     QUANTILE TO BE ESTIMATED
C           ALPHA     =     PARAMETER OF G1 DISTRIBUTION
C
C     N.B.  DDGAM IS CONTAINED IN VARMOM SUBROUTINE FILE
C
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION PQ,Q1P,ALPHA
      
         Q1P      =   FP_G1_ICDF(PQ,ALPHA)
         D_G1_INV =  -DDGAM(ALPHA,Q1P)/FP_G1_PDF(Q1P,ALPHA)
         
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
         SUBROUTINE EXPMOMDERIV(PARMS,XMIN,XMAX,JAC)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE DERIVATIVE OF EXPECTED VALUE OF 
C         MOMENTS OF GAMMA VARIATE TRUNCATED BETWEEN XMIN AND XMAX
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 10, 1999
C
C     REQUIRES FUNCTIONS FOR LOG-GAMMA AND PSI
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION 
     1  PARMS(3),XMIN,XMAX,JAC(3,3),
     2  JAC1(3,3),JAC2(3,3),ALPHA,BETA,TAU,M(3)

          TAU   = PARMS(1)
          ALPHA = PARMS(2)
          BETA  = PARMS(3)
        CALL DEXPECT(TAU,ALPHA,BETA,XMIN,XMAX,M,JAC1)
         
        CALL DPDM(PARMS,JAC2)
        
        CALL DMRRRR(3,3,JAC1,3,3,3,JAC2,3,3,3,JAC,3)
        
        RETURN
        END
        
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION DDGAM(ALPHA,X)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE DERIVATIVE OF INCOMPLETE GAMMA FUNCTION
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 9, 1999
C
C     REQUIRES FUNCTIONS FOR LOG-GAMMA AND PSI
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION I,TOL,LOGX,SUM
      
       DATA TOL/1.D-11/
      
       IF(ALPHA .LE. 0.D0 .OR. X .LE. 0.D0) THEN
               RESULT = 0.D0
               GOTO 100
       END IF
              
       IF(ABS(X-ALPHA)/SQRT(ALPHA+1.D0) .GT. 7.D0) THEN
               RESULT = 0.D0
               GOTO 100
       END IF
      
          A        = ALPHA
          LOGX     =  LOG(X)
C
      IF(X .LT. 5.D0) THEN
          T        =  X**A
          R        =  (A*LOGX-1.D0)/A**2
          SUM      =  T*R
        DO 10 I2=1,1000
          I        =  I2
          AI       =  A+I
          T        =  -T*X/I
          R        =  (AI*LOGX-1.D0)/AI**2
          DEL      =  R*T
          SUM      =  SUM + DEL
           IF(I .GT. 1 .AND. ABS(DEL) .LT. (1.D0+ABS(SUM))*TOL) THEN
             RESULT =  SUM/EXP(DGAMLN(A,IERR)) - DPSI(A)*FP_G1_CDF(X,A)
             GOTO 100
           ENDIF
10      CONTINUE
          GOTO 99
C
      ELSE IF(X .GT. A+30.D0) THEN
          T        =  EXP(-X + (A-1.D0)*LOGX - DGAMLN(A,IERR) )
          R        =  LOGX - DPSI(A)
          SUM      =  R*T
        DO 20 I2=1,INT(A-1.D0)
          I         =  I2
          AMI       =  A-I
          T        =  T*(AMI)/X
          R        =  LOGX - DPSI(AMI)
          DEL      =  R*T
          SUM      =  SUM + DEL
           IF(I .GT. 1 .AND. ABS(DEL) .LT. (1.D0+ABS(SUM))*TOL) THEN
             RESULT = -SUM
             GOTO 100
           ENDIF
20      CONTINUE
          GOTO 99
C
      ELSE
          T        =  EXP(-X + A*LOGX - DGAMLN(A+1.D0,IERR) )
          R        =  LOGX - DPSI(A+1.D0)
          SUM      =  R*T
        DO 30 I2=1,10000
          I        =  I2
          T        =  T*X/(A+I)
          R        =  LOGX - DPSI(A+I+1.D0)
          DEL      =  R*T
          SUM      =  SUM + DEL
           IF(I .GT. 1 .AND. ABS(DEL) .LT. (1.D0+ABS(SUM))*TOL) THEN
             RESULT = SUM
             GOTO 100
           ENDIF
30      CONTINUE
          GOTO 99
      ENDIF
C
100   CONTINUE
        DDGAM      =  RESULT
      RETURN
99    CONTINUE
              WRITE(*,*) 'SOMETHING BAD (DDGAM).....'
              WRITE(*,*) 'Contact Tim Cohn (703/648-5711).....'
              WRITE(*,*) 'tacohn@usgs.gov.....'
              WRITE(*,*) 'ALPHA: ',ALPHA
              WRITE(*,*) 'X:     ',X
              READ(*,*)
              STOP
      END

C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      SUBROUTINE DEXPECT(TAU,ALPHA,BETA,TL,TU,M3,DM3)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE DERIVATIVE OF INCOMPLETE GAMMA FUNCTION
C       TRUNCATED BELOW AT TL AND ABOVE AT TU, WITH PARMS TAU, ALPHA, BETA
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 9, 1999
C
C     REQUIRES FUNCTIONS FOR LOG-GAMMA AND PSI
C
C     NOTE RETURN VARIABLES:
C         M3[J]      CONTAINS E[X^J]
C         DM3[J,K]   CONTAINS D[ E[X^J] ]/P[K]
C
C          DM3[1,1]  =   D[ E[X] ]/D_TAU
C          DM3[1,2]  =   D[ E[X] ]/D_ALPHA
C          DM3[1,3]  =   D[ E[X] ]/D_BETA
C
C          DM3[2,1]  =   D[ E[X^2] ]/D_TAU
C          DM3[2,2]  =   D[ E[X^2] ]/D_ALPHA
C          DM3[2,3]  =   D[ E[X^2] ]/D_BETA
C
C          DM3[3,1]  =   D[ E[X^3] ]/D_TAU
C          DM3[3,2]  =   D[ E[X^3] ]/D_ALPHA
C          DM3[3,3]  =   D[ E[X^3] ]/D_BETA
C
C       WHERE P[1] = TAU, P[2] = ALPHA, P[3] = BETA
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      SAVE
      DOUBLE PRECISION 
     1    ALPHA,BETA,TAU,TL,TU,
     2    DG1(0:3),G1(0:3),
     3    ADJ(0:3),DADJ(0:3),B(0:3),DB(0:3),
     4    M1(3),DM1(0:3,3),M2(3),DM2(3,3),M3(3),DM3(3,3),
     5    LSTDL,LSTDU,FU(0:3),FL(0:3)

C
C    RETURN IF THE RESULTS DO NOT DEPEND ON THE PARAMETERS
C
        IF(TU .LE. TL) GOTO 99
        
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C    STANDARDIZE THE THRESHOLDS;  COMPUTE DERIVATIVES
C
        IF(BETA .GT. 0.D0) THEN
          STDL  =  MAX(0.D0,(TL-TAU)/BETA)
          STDU  =  (TU-TAU)/BETA
        ELSE IF(BETA .LT. 0.D0) THEN
          STDL  =  MAX(0.D0,(TU-TAU)/BETA)
          STDU  =  (TL-TAU)/BETA
        ELSE IF(BETA .EQ. 0.D0) THEN
          GOTO 99
        ENDIF
          LSTDL   =  LOG(MAX(1.D-99,STDL))
          LSTDU   =  LOG(MAX(1.D-99,STDU))
          
          DSTDLT  =  -1.D0/BETA
          DSTDUT  =  -1.D0/BETA
          DSTDLB  =  -STDL/BETA
          DSTDUB  =  -STDU/BETA

C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C    1.  FIRST FOR THE ONE-PARAMETER GAMMA (ASSUME TAU=0; BETA=1)
C         -- NOTE THE FOLLOWING:
C              E[X^P] = GAMMA[ALPHA+P,TL,TU] / GAMMA[ALPHA,TL,TU]
C                     = ADJ[P]*B[P]
C              ADJ[P] = GAMMA[ALPHA+P]/GAMMA[ALPHA]
C              B[P]   = GCDF[ALPHA+P,TL,TU] / GCDF[ALPHA,TL,TU]
C              GCDF[ALPHA,TL,TU] = GAMMA[ALPHA,TL,TU]/GAMMA[ALPHA,0,INF]
C

C
C      1.A   FIRST COMPUTE THE ADJUSTMENT FACTOR
C
          ADJ(0)    =  1.D0
          DADJ(0)   =  0.D0
        DO 40 J=1,3
          ADJ(J)    =  (ALPHA+J-1.D0) * ADJ(J-1)
          DADJ(J)   =  ADJ(J-1) + (ALPHA+J-1.D0) * DADJ(J-1)
40      CONTINUE

C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C      1.B   SECOND COMPUTE THE INCOMPLETE G1 CDF AND ITS DERIVATIVES
C              -- NB:  B(J) IS THE NORMALIZED CDF
C                      DB(J) IS COMPUTED FROM RATIO
C                      G1 CONTAINS TAU, BETA THROUGH STDL
C
      DO 20 J=0,3
        GU      =  FP_G1_CDF(STDU,ALPHA+J)
        GL      =  FP_G1_CDF(STDL,ALPHA+J)
        G1(J)   =  GU-GL
        DG1(J)  =  DDGAM(ALPHA+J,STDU) - DDGAM(ALPHA+J,STDL)
        B(J)    =  G1(J)/G1(0)
        DB(J)   =  (DG1(J)*G1(0) - DG1(0)*G1(J))/G1(0)**2
        FU(J)   =  FP_G1_PDF(STDU,ALPHA+J)
        FL(J)   =  FP_G1_PDF(STDL,ALPHA+J)
        DM1(J,1) = (
     1     (FU(J)*DSTDUT - FL(J)*DSTDLT) * G1(0) -
     2     (FU(0)*DSTDUT - FL(0)*DSTDLT) * G1(J) )/G1(0)**2
        DM1(J,3) = (
     1     (FU(J)*DSTDUB - FL(J)*DSTDLB) * G1(0) -
     2     (FU(0)*DSTDUB - FL(0)*DSTDLB) * G1(J) )/G1(0)**2
20    CONTINUE
    
C
C      1.C   NOW COMPUTE THE DERIVATIVE OF THE PRODUCT
C
      DO 50 J=1,3
        M1(J)     =  ADJ(J)*B(J)
        DM1(J,2)  =  ADJ(J)*DB(J) + DADJ(J)*B(J)
        DM1(J,1)  =  ADJ(J)*DM1(J,1)
        DM1(J,3)  =  ADJ(J)*DM1(J,3)
50    CONTINUE

C
C
C    2.  NOW GENERALIZE TO THE 2-PARAMETER GAMMA (ASSUME TAU=0)
C        --NOTE THAT ORDER OF PARAMETERS IS TAU,ALPHA,BETA
C          (TAC -- ANOTHER GOOFY CONVENTION); THEREFORE
C              DM2(J,1) REFERS TO TAU
C              DM2(J,2) REFERS TO ALPHA
C              DM2(J,3) REFERS TO BETA
C
      DO 60 J=1,3
        M2(J)     = BETA**J * M1(J)
        DM2(J,1)  = BETA**J * DM1(J,1)
        DM2(J,2)  = BETA**J * DM1(J,2)
        DM2(J,3)  = J * BETA**(J-1) * M1(J) + BETA**J * DM1(J,3)
60    CONTINUE
C
C
C    3.  NOW GENERALIZE TO THE 3-PARAMETER GAMMA
C
C
      DO 70 J=1,3
          M3(J)     = M2(J) + TAU**J
          DM3(J,1)  = J * POWER(TAU,(J-1)) + DM2(J,1)
          DM3(J,2)  = DM2(J,2)
          DM3(J,3)  = DM2(J,3)
        DO 70 K=1,J-1
          CH        =  CHOOSE(J,K)
          M3(J)     =  M3(J) + CH*TAU**(J-K) * M2(K)
          DM3(J,1)  =  DM3(J,1) + CH*(
     1                   (J-K)*POWER(TAU,(J-K-1)) * M2(K) +
     2                   TAU**(J-K) * DM2(K,1) )
          DM3(J,2)  =  DM3(J,2) + CH*TAU**(J-K) * DM2(K,2)
          DM3(J,3)  =  DM3(J,3) + CH*TAU**(J-K) * DM2(K,3)
70    CONTINUE

      RETURN
C 
C  NON-FUNCTIONAL EXIT
C
99    CONTINUE
          CALL DSET(9,0.D0,DM3,1)
          CALL DSET(3,0.D0,M3,1)
          RETURN
      END
C
      DOUBLE PRECISION FUNCTION POWER(X,I)
C      SAVE
      DOUBLE PRECISION X
      INTEGER I
      IF(X .EQ. 0.D0 .AND. I .EQ. 0) THEN
         POWER = 1.D0
      ELSE
         POWER = X**I
      ENDIF
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_TNC_CDF(TP,NU,DELTA)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE NON-CENTRAL T CDF
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 9, 1999
C
C     REQUIRES FUNCTIONS FOR Z
C
C     REFERENCE:  APPROXIMATE FORMULA ON PAGE 949, A&S
C                 TNC FUNCTION IS ALGORITHM AS 243  
C                    APPL. STATIST. (1989), VOL.38, NO. 1
C
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT NONE
C      SAVE
      
      INTEGER IFAULT
      
      DOUBLE PRECISION TP,NU,DELTA,Z,FP_Z_CDF,TNC,ANS
      
      IF( NU .GT. 20.D0) THEN
       Z  = (TP*(1.D0-1.D0/(4.D0*NU))-DELTA)/SQRT(1.D0+TP**2/(2.D0*NU))
       ANS=  FP_Z_CDF(Z) 
      ELSE
       ANS=  TNC(TP,NU,DELTA,IFAULT)
       IF(IFAULT .NE. 0) THEN
          WRITE(*,*) 'TROUBLE IN FP_TNC_CDF:TP,NU,DELTA',TP,NU,DELTA
       Z  = (TP*(1.D0-1.D0/(4.D0*NU))-DELTA)/SQRT(1.D0+TP**2/(2.D0*NU))
       ANS=  FP_Z_CDF(Z) 
       ENDIF
      ENDIF
      
      FP_TNC_CDF =  ANS
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_TNC_ICDF_OLD(P,NU,DELTA)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE NON-CENTRAL T INVERSE CDF
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 9, 1999
C
C     REQUIRES FUNCTIONS FOR Z-INVERSE
C
C     REFERENCE:  INVERSION OF APPROXIMATE FORMULA ON PAGE 949, A&S
C
C     NOTE:  THIS IS NOT VERY ACCURATE (TAC: 20 SEP 2010)
C            USED TO BE CALLED FP_TNC_ICDF; NOW FP_TNC_ICDF_OLD
C            FORMER FP_TNC2_ICDF NOW FP_TNC_ICDF (TAC: 20 SEP 2010)
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT NONE
C      SAVE
      DOUBLE PRECISION P,NU,DELTA,X,ANS,FP_Z_ICDF,FP_TNC_ICDF

      IF(NU .LE. 20.D0) THEN
           ANS = FP_TNC_ICDF(P,NU,DELTA)
      ELSE
      X    =  FP_Z_ICDF(P)
      ANS  =  (4.D0*NU*(DELTA*(-1.D0 + 4.D0*NU) + 
     1           X*SQRT(1.D0 + 16.D0*NU**2 + 
     2           8.D0*NU*(-1.D0 + DELTA**2 - X**2))))/
     3           (1.D0 + 16.D0*NU**2 - 8.D0*NU*(1.D0 + X**2))
      ENDIF
     
      FP_TNC_ICDF_OLD =  ANS
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FP_TNC_ICDF(P_IN,NU_IN,DELTA_IN)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE NON-CENTRAL T INVERSE CDF
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 9, 1999
C
C     REQUIRES FUNCTIONS FOR Z-INVERSE, GAMMA INVERSE
C
C     NOTE:  THIS USED TO BE LABELED FP_TNC2_ICDF (TAC: 20 SEP 2010)
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT NONE
C      SAVE
      
      INTEGER IFLAG
      
      DOUBLE PRECISION 
     1   P_IN,NU_IN,DELTA_IN,
     2   P,NU,DELTA,
     3   AE,RE,B,C
     
      DOUBLE PRECISION FTNCINV
      
      EXTERNAL FTNCINV
      
      COMMON /ZNCT02/P,NU,DELTA

      DATA RE/1.D-8/,AE/1.D-8/

      P     =  P_IN
      NU    =  NU_IN
      DELTA =  DELTA_IN
      
c  change made to limits 20 Sep 2010 (TAC)
c      X1    =  FP_Z_ICDF(P)+DELTA
c                PARMS(1)  =  NU/2.D0
c                PARMS(2)  =  2.D0
c      X2    =  SQRT(FP_G2_ICDF(1.D0-P,PARMS)/NU)
c      
c      IF(X2 .LE. 0.D0) GOTO 99
c      
c        B     =  MIN(0.D0,X1,1.D0/X2,X1/X2)
c        C     =  MAX(0.D0,X1,1.D0/X2,X1/X2)

         B     =  -32.D0  !  Works for 0.01<=p<=0.99,nu=1,nc=0  07/15/2013 (TAC)
         C     =   32.D0 
                      
      CALL DZEROIN(FTNCINV,B,C,RE,AE,IFLAG)
      
      IF(IFLAG .NE. 1) GOTO 99
      
      FP_TNC_ICDF =  B
      
      RETURN
      
99    CONTINUE
        WRITE(*,*) 'ERROR IN FP_TNC_ICDF',IFLAG,P,NU,DELTA
C        READ(*,*)
        FP_TNC_ICDF =  0.D0
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
      DOUBLE PRECISION FUNCTION FTNCINV(TP)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE CDF OF NON-CENTRAL T DISTRIBUTION
C
C     AUTHOR.....TIM COHN
C     DATE.......APRIL 19, 1999
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c        SAVE
        DOUBLE PRECISION TP,ANS,P,NU,DELTA,FP_TNC_CDF
      
        COMMON /ZNCT02/P,NU,DELTA

          ANS     =    FP_TNC_CDF(TP,NU,DELTA) - P
          FTNCINV =    ANS
          
        RETURN
        END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C
      DOUBLE PRECISION FUNCTION FP_B_CDF(X,A,B)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF BETA(X,ALPHA,BETA)
C
C     AUTHOR.....TIM COHN
C     DATE.......08 MAR 2010
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      DOUBLE PRECISION
     1  P,Q,X,A,B,BOUND
      
      INTEGER
     1  WHICH,STATUS
       
      DATA WHICH/1/

      CALL CDFBET ( WHICH, P, Q, X, 1.D0-X, A, B, STATUS, BOUND )
        IF(STATUS .NE. 0) WRITE(*,*) ' ERROR IN CDFBET/FP_B_CDF',STATUS
      
      FP_B_CDF = P
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C
      DOUBLE PRECISION FUNCTION FP_B_ICDF(P,A,B)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE INVERSE CDF OF BETA(P,A,B)
C
C     AUTHOR.....TIM COHN
C     DATE.......08 MAR 2010
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      DOUBLE PRECISION
     1  P,X,Y,A,B,BOUND
      
      INTEGER
     1  WHICH,STATUS
       
      DATA WHICH/2/

      CALL CDFBET ( WHICH, P, 1.D0-P, X, Y, A, B, STATUS, BOUND )
        IF(STATUS .NE. 0) WRITE(*,*) ' ERROR IN CDFBET/FP_B_ICDF',STATUS
      
      FP_B_ICDF = X
      
      RETURN
      END
c
c
c  Note:  The following two commented-out functions do not work
c         They do not deal with the non-centrality parameter at all
c
c         tac 15 Sep 2010
c
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C
      DOUBLE PRECISION FUNCTION FP_T_CDF(T,DF)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF CENTRAL T(T,DF)
C
C     AUTHOR.....TIM COHN
C     DATE.......08 MAR 2010
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      DOUBLE PRECISION
     1  DF,P,Q,T,BOUND
      
      INTEGER
     1  WHICH,STATUS
       
      DATA WHICH/1/

      CALL CDFT ( WHICH, P, Q, T, DF, STATUS, BOUND )
        IF(STATUS .NE. 0) WRITE(*,*) ' ERROR IN CDFT/FP_T_CDF2',STATUS
              
      FP_T_CDF = P
      
      RETURN
      END
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
C
      DOUBLE PRECISION FUNCTION FP_T_ICDF(P,DF)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF CENTRAL T(T,DF)
C
C     AUTHOR.....TIM COHN
C     DATE.......08 MAR 2010
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      DOUBLE PRECISION
     1  DF,P,T,BOUND
      
      INTEGER
     1  WHICH,STATUS
       
      DATA WHICH/2/
 
        CALL CDFT ( WHICH, P, 1.D0-P, T, DF, STATUS, BOUND )
        IF(STATUS .NE. 0) WRITE(*,*) ' ERROR IN CDFT/FP_T_CDF2',STATUS
      
      FP_T_ICDF = T
      
      RETURN
      END
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      INTEGER FUNCTION MGBTP(X,N,PVALUEW)
!*** REVISION 2012.09
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
! 
!       IDENTIFIES THE NUMBER OF LOW OUTLIERS USING THE 
!       GENERALIZED GRUBBS-BECK TEST
!       NOTES:  SEE 2011 MANUSCRIPT BY COHN, STEDINGER, ENGLAND, ET AL.
! 
!       AUTHOR.........TIM COHN
!       DATE...........1 MAR 2010 (TAC)
!               MODIFIED....... 8 MAR 2010 (TAC)
!               MODIFIED.......29 APR 2010 (TAC) [OFFSET ADDED FOR SMALL N]
!               MODIFIED.......21 OCT 2010 (TAC) [OFFSET REMOVED; ESTIMATOR
!                                  FOR E[S] REPLACED WITH GAMMA FUNCTION
!               MODIFIED.......5 JAN 2011 (TAC)
!               MODIFIED.......16 FEB 2011 (TAC) RETURNS PVALUEW
!                 MODIFICATIONS FOR MGBT TESTING
!                  John England, Bureau of Reclamation, 17-SEP-2012                     
!               MODIFIED.......26 SEP 2012 (TAC) REWORKED CODE TO CLARIFY LOGIC
!                                                AND AVOID STDOUT CALLS
!
!c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
!       calling routine: gbtest
!
!       input variables (from subroutine gbtest):
!       ------------------------------------------------------------------------
!            n          i*4  number of observations (uncensored, with ql=qu)
!            x(n)       r*8  vector of floods in real space (qs from gbtest)
!
!                       note: see subroutine gbtest for input and call
!                       input is systematic 'Syst' records where ql_in=qu_in
!                       and tl_in <= gage_base from peakfqSA
!
!****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
!       output variables:
!       ------------------------------------------------------------------------
!           pvaluew(n)  r*8  vector of MGB p-values
!                           (function of n,k,w) returned from ggbcritp
!
!****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
!       controlling parameters:
!       ------------------------------------------------------------------------
!           Alphaout    r*8  significance level for outward sweep from median
!                            valid value range: >0
!
!           Alphain     r*8  significance level for inward sweep from smallest obs
!                            valid value range: any; controlled by tiny
!
!           Alphazeroin r*8  significance level for inward sweep from smallest obs
!                            valid value range: >0
!
!       default parameter values stored in block data
!         data Alphaout/0.005d0/,Alphain/0.d0/,Alphazeroin/0.1d0/      
!
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
      IMPLICIT NONE
      
      INTEGER DIM1,N
C JFE      PARAMETER (DIM1=10000)
      PARAMETER (DIM1=20000)
      
      DOUBLE PRECISION 
     * X(N),PVALUEW(DIM1),W(DIM1),ZT(DIM1),S1,S2,XM,XV,
     1 GGBCRITP,Alphaout,Alphain,Alphazeroin

      INTEGER 
     * I,J1,J2,J3,N2,NC

      common /MGB001/Alphaout,Alphain,Alphazeroin

      DO 10 I=1,N
         ZT(I) = LOG10(MAX(1D-88,X(I)))
         PVALUEW(I) = -99.D0    ! MISSING VALUES
 10   CONTINUE

!     sort log flows from smallest to largest
      CALL DSVRGN(N,ZT,ZT)
        
!     set starting point for MGBT search at approximate median position (1/2 N)
      N2    = N/2

      S1 = 0.D0
      S2 = 0.D0

      DO 20 I=N,N2+2,-1
         S1 = S1 + ZT(I)
         S2 = S2 + ZT(I)**2
 20   CONTINUE

      DO 30 I=N2,1,-1
         S1  = S1 + ZT(I+1)
         S2  = S2 + ZT(I+1)**2
         NC  = N-I
         XM  = S1/DBLE(NC)
         XV  = (S2-NC*XM**2)/(NC-1.D0)
         W(I)       = (ZT(I)-XM)/SQRT(XV)
         PVALUEW(I) = GGBCRITP(N,I,W(I))
 30   CONTINUE
!
!     Determine number of low outliers in 2 or 3 steps.
!     Based on TAC original code and JRS recommendations.
!
!     Step 1. Outward sweep from median (always done).
!             alpha level of test = Alphaout
!             number of outliers = J1
!
!     Step 2. Inward sweep from largest low outlier identified in Step 1.
!             alpha level of test = Alphain
!             number of outliers = J2
!
!     Step 3. Inward sweep from smallest observation
!             alpha level of test = Alphazeroin
!             number of outliers = J3
!
!     Initialize counters
      J1 = 0  ! Outward sweep number of low outliers
      J2 = 0  ! Inward sweep number of low outliers
      J3 = 0  ! Oth Inward sweep number of low outliers

!     1) Outward sweep check: Loop over low flows up to median
      DO 110 I=N2,1,-1
         IF(PVALUEW(I) .LT. Alphaout) GOTO 111
110   CONTINUE
111   CONTINUE
        J1=I

!     2) Inward sweep check with Alphain
           J2 = J1
      DO 120 I=J1+1,N2  
        IF( (PVALUEW(I) .GE. Alphain)) GOTO 121
120   CONTINUE
121   CONTINUE
        J2 = I-1

!     3) Inward sweep check with Alphazeroin
      DO 130 I=1,N2
        IF( (PVALUEW(I) .GE. Alphazeroin)) GOTO 131
130   CONTINUE
131   CONTINUE
        J3=I-1
        
!     Set  number of low outliers as max of 3 sweeps
      MGBTP = MAX(J1,J2,J3)
      RETURN
      END
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      blockdata MGBblk
      implicit none
      double precision Alphaout,Alphain,Alphazeroin     
      common /MGB001/Alphaout,Alphain,Alphazeroin      
      data Alphaout/0.005d0/,Alphain/0.d0/,Alphazeroin/0.1d0/      
      end
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      INTEGER FUNCTION MGBT(X,N)
!*** REVISION 2011.01
!  RETURNS JUST NUMBER OF LOW OUTLIERS; MGBTP ALSO RETURNS PVALUES
      IMPLICIT NONE
      INTEGER DIM1,N,MGBTP
C JFE      PARAMETER (DIM1=10000)
      PARAMETER (DIM1=20000)
      DOUBLE PRECISION X(*),PVALUEW(DIM1)
      MGBT = MGBTP(X,N,PVALUEW)
      RETURN
      END

        DOUBLE PRECISION FUNCTION GGBCRITP(N,R,ETA)
!*** REVISION 2011.01
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
! 
!       PROGRAM TO COMPUTE P-VALUES (GGCRITP) AND CRITICAL POINTS  (GGCRITK)
!         FOR A GENERALIZED GRUBBS-BECK TEST
!       NOTES:  SEE 2010 MANUSCRIPT BY COHN, ENGLAND, ET AL.
! 
!       AUTHOR.........TIM COHN
!       DATE...........1 MAR 2010 (TAC)
!               MODIFIED....... 8 MAR 2010 (TAC)
!               MODIFIED.......29 APR 2010 (TAC) [OFFSET ADDED FOR SMALL N]
!               MODIFIED.......21 OCT 2010 (TAC) [OFFSET REMOVED; ESTIMATOR
!                                  FOR E[S] REPLACED WITH GAMMA FUNCTION
!               MODIFIED........5 JAN 2011 (TAC) MGBT REVISED
!
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
      IMPLICIT NONE
      
      INTEGER
     1 LIMIT
     
      PARAMETER
     1 (LIMIT=40000)
     
      INTEGER
     * N,R,N_IN,R_IN,
     1 NEVAL,IER,LENW,LAST,IWORK(LIMIT),
     2 KEY
     
      DOUBLE PRECISION 
     * ETA,ETA_IN,
     1 FGGB,EPSABS,EPSREL,RESULT,ABSERR,WORK(4*LIMIT)
     

      EXTERNAL FGGB
      
      COMMON /GGB001/ETA_IN,N_IN,R_IN
 
      DATA EPSABS/1.D-5/,EPSREL/1.D-5/,KEY/3/,LENW/160000/
      
     
      IF(N .LT. 10 .OR. R .GT. N/2) THEN
        WRITE(*,*) ' INVALID INPUT DATA/GGBCRITP (N,R): ', N,R
        GGBCRITP = 0.5D0
        RETURN
      ELSE
          N_IN   = N
          R_IN   = R
          ETA_IN = ETA
      ENDIF
      
      CALL DQAG(FGGB,0.D0,1.D0,EPSABS,EPSREL,KEY,RESULT,ABSERR,
     *    NEVAL,IER,LIMIT,LENW,LAST,IWORK,WORK)

      GGBCRITP   =  RESULT
      
      RETURN
      END
            
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      DOUBLE PRECISION FUNCTION FGGB(PZR)
!*** REVISION 2011.01
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
      IMPLICIT NONE
      
      INTEGER
     * N,R
     
      DOUBLE PRECISION 
     * PZR,ETA,DF,MuM, MuS2,VarM,VarS2,CovMS2,EX1,EX2,EX3,EX4,
     2 CovMS,VarS,alpha,beta,dlngam,
     3 MuMP,EtaP,H,Lambda,MuS,ncp,q,VarMP,PR,ZR,N2,
     4 FP_B_ICDF,FP_Z_ICDF,FP_Z_PDF,
     5 FP_TNC_CDF

      COMMON /GGB001/ETA,N,R
!
!  Compute the value of the r-th smallest obs. based on its order statistic
!
      N2   = dble(N-R)
      PR   = FP_B_ICDF(PZR,dble(R),dble(N+1-R))
      ZR   = FP_Z_ICDF(PR)

!
!  Calculate the expected values of M, S2, S and their variances/covariances
!
      H        =  FP_Z_PDF(ZR)/(Max(1.d-10,1.d0-PR))
      EX1      =  H
      EX2      =  1 + H * ZR
      EX3      =  2 * EX1 + H * ZR**2
      EX4      =  3 * EX2 + H * ZR**3
      
      MuM      =  EX1
      MuS2     =  (EX2-EX1**2)
      VarM     =  MuS2/N2
      VarS2    =  ((EX4 - 4*EX3*EX1 + 6*EX2*EX1**2 - 3*EX1**4) - 
     1                MuS2**2)/N2 + 2.d0/((N2-1.d0)*n2)*MuS2**2

       alpha    =  MuS2**2/VarS2
       beta     =  MuS2/alpha

      CovMS2   =  (EX3 - 3*EX2*EX1 + 2*EX1**3)/SQRT(N2*(N2-1.D0))
      
      MuS      =  SQRT(beta) * exp(dlngam(alpha+0.5d0)-dlngam(alpha))
      CovMS    =  CovMS2/(2*MuS)
      VarS     =  MuS2 - MuS**2
      
        Lambda = CovMS/VarS
        EtaP   = Eta + Lambda
        MuMP   = MuM - Lambda * MuS
        VarMP  = VarM - CovMS**2/VarS
        df     = 2.d0*alpha
        ncp    = (MuMP -zr)/SQRT(VarMP)
        q      = -sqrt(MuS2/VarMP) * EtaP
       FGGB = 1.d0 - FP_TNC_CDF(q,df,ncp)
c        FGGB = 1.d0 - FP_NonCentral_T_CDF(q,df,ncp)

      return
      end
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!*** 
!*** 
      DOUBLE PRECISION FUNCTION GGBCRITK(N,R,P)
!*** REVISION 2011.01
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
! 
!       PROGRAM TO COMPUTE CRITICAL POINTS FOR A GENERALIZED GRUBBS-BECK
!         TEST
!       NOTES:  SEE 2010 MANUSCRIPT BY COHN, ENGLAND, ET AL.
! 
!       AUTHOR.........TIM COHN
!       DATE...........25 OCT 2010 (TAC)
!               MODIFIED....... 8 MAR 2010 (TAC)
!
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
      IMPLICIT NONE
      
CTAC      INTEGER
CTAC     1 LIMIT
CTAC     
CTAC      PARAMETER
CTAC     1 (LIMIT=40000)
     
      INTEGER
     * N,R,N_IN,R_IN,IFLAG,NCALLS
     
      DOUBLE PRECISION 
     * FGGBX,FP_Z_ICDF,
     1 P,P_IN,B,C,GUESS,RE,AE

      EXTERNAL FGGBX
      
      COMMON /GGB002/P_IN,N_IN,R_IN,NCALLS
 
      DATA RE/1.D-6/,AE/1.D-6/
      
      IF(R .LT. 1 .OR. R .GT. (N-1)) THEN
        WRITE(*,*) 'ERROR IN GGBCRITK, R,N: ',R,N
        STOP
      ENDIF
       N_IN  = N
       R_IN  = R
       P_IN  = P
       GUESS = FP_Z_ICDF(DBLE(R)/DBLE(N))
         NCALLS = 0
         B = -20.d0
         C =  20.d0
      CALL DFZERO (FGGBX, B, C, GUESS, RE, AE, IFLAG)
        IF(IFLAG .NE. 1) THEN
          WRITE(*,*) 'ERROR DFZERO/GGBCRITK, IFLAG: ',IFLAG
          WRITE(*,*)
        ENDIF
        
      GGBCRITK = B
      RETURN
      END
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      DOUBLE PRECISION FUNCTION FGGBX(ETA)
!*** REVISION 2011.01
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
      IMPLICIT NONE
      
      INTEGER
     * N,R,NCALLS
     
      DOUBLE PRECISION 
     * GGBCRITP,ETA,P

      COMMON /GGB002/P,N,R,NCALLS
!
!  Compute the value of the r-th smallest obs. based on its order statistic
!
      FGGBX = GGBCRITP(N,R,ETA) - P
        NCALLS = NCALLS+1
      RETURN
      END
      
      DOUBLE PRECISION FUNCTION FP_NonCentral_T_CDF(X,DF_IN,NCP_IN)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C 
C     PROGRAM TO COMPUTE CDF OF NONCENTRAL T(T,DF,NCP)
C
C     AUTHOR.....TIM COHN
C     DATE.......01 APR 2010
C
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////

      IMPLICIT NONE
      
      INTEGER
     1 LIMIT
     
      PARAMETER
     1 (LIMIT=40000)
     
      INTEGER
     1 NEVAL,IER,LENW,LAST,IWORK(LIMIT),KEY
CTAC
CTAC

      DOUBLE PRECISION 
     1 BOUND,EPSABS,EPSREL,RESULT,ABSERR,WORK(4*LIMIT),
     2 NCP,W

CTAC
      DOUBLE PRECISION
     1  DF_IN,NCP_IN,X,DF
      
      INTEGER
     1  INF
     
      COMMON /FPNC01/W,DF,NCP

      DATA BOUND/0.D0/,INF/1/
      
      DATA EPSABS/1.D-10/,EPSREL/1.D-10/,KEY/3/,LENW/160000/
     
      EXTERNAL FNCT01
      
      W   =  X
      DF  =  DF_IN
      NCP =  NCP_IN
       
      CALL DQAGI(FNCT01,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,
     1   IER,LIMIT,LENW,LAST,IWORK,WORK)
      
      FP_NonCentral_T_CDF = RESULT
      
      RETURN
      END
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      DOUBLE PRECISION FUNCTION FNCT01(S2)
!*** REVISION 2011.01
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
!****|==|====-====|====-====|====-====|====-====|====-====|====-====|==////////
!
      IMPLICIT NONE
      
      DOUBLE PRECISION 
     * W,DF,NCP,S2,PARMS(2),
     2 FP_Z_CDF,FP_G2_PDF

      COMMON /FPNC01/W,DF,NCP
      
        PARMS(1) = DF/2.D0
        PARMS(2) = 2.D0/DF

      FNCT01 = FP_Z_CDF(W*SQRT(S2)-NCP) * FP_G2_PDF(S2,PARMS)
      
      RETURN
      END 
c*_*-*~*_*-*~*             new program begins here            *_*-*~*_*-*~*_*-*~
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c****|subroutine plotposHS
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c    
c    this routine provides plotting positions for censored data 
c      using the Hirsch/Stedinger plotting position formula (hirsch, 1987)
c      Purpose is operational version of peakfq (bulletin 17b
c      implementation) which would include the expected moments
c      algorithm (ema; cohn et al. 1997; 2001)
c
c    N.B:  The returned values are "exceedance probabilities" which are 
c             equal to 1-NonExceedance probabilities
c       
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c    
c    References
c    ----------
c
c    @article{hirsch1987plotting,
c      title={{Plotting positions for historical floods and their precision}},
c      author={Hirsch, R.M. and Stedinger, J.R.},
c      journal={Water Resources Research},
c      volume={23},
c      number={4},
c      pages={715--727},
c      year={1987}
c    }
c
c    @article{cohn1997algorithm,
c      title={{An algorithm for computing moments-based flood quantile estimates when historical flood information is available}},
c      author={Cohn, TA and Lane, WL and Baier, WG},
c      journal={Water Resources Research},
c      volume={33},
c      number={9},
c      pages={2089--2096},
c      year={1997}
c    }
c
c    @article{cohn2001confidence,
c      title={{Confidence intervals for Expected Moments Algorithm flood quantile estimates}},
c      author={Cohn, T.A. and Lane, W.L. and Stedinger, J.R.},
c      journal={Water resources research},
c      volume={37},
c      number={6},
c      pages={1695--1706},
c      year={2001}
c    }
c
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    development history
c
c    timothy a. cohn        15 apr 2010
c       modified            16 apr 2010
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    subroutine calls
c
c    ppfit                  available in probfun.f
c    arrange                available in probfun.f
c    porder                 available in probfun.f
c    ppsolve                available in probfun.f
c    plpos                  available in probfun.f
c    shellsort              available in probfun.f
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       n.b.  use of plotposHS requires two distinct types of information:
c
c            a.  "observations," which are described by {ql,qu};
c                  specifically, the true value of q is assumed to lie
c                  wihin the interval ql(i) <= q(i) <= qu(i);
c                  if q(i) is known exactly, then ql(i) = q(i) = qu(i);
c
c                  {ql,qu} are required by ema to fit p3 distn to data
c
c            b.  "censoring pattern," which is described by {tl,tu}
c                  {tl,tu} define the interval on which data get observed:
c                  if q is not inside {tl,tu}, q is either left or right 
c                  censored.
c
c                  examples:
c                    {-inf,+inf}  => systematic data (exact value of q
c                      would be measured)
c                    {t,   +inf}  => historical data 
c                      q>=t would be known exactly; 
c                      q<t would be reported as censored ("less-than")
c                
c                    {tl,tu} are required by ema to estimate confidence 
c                      intervals; they are not required for frequency
c                      curve determination
c
c       n.b.  broken record is represented in the time-series data by a
c             threshold value, tl(i), equal to 1.d12 or larger.  This
c             indicates that only flows greater than 1.d12 would have 
c             been recorded.  No such floods have ever occurred on the
c             planet [TAC, 2010, personal communication].
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       input variables:
c       ------------------------------------------------------------------------
c            n          i*4  number of observations (censored, uncensored, or 
c                              other)
c            ql(n)      r*8  vector of lower bounds on floods
c            qu(n)      r*8  vector of upper bounds on floods
c            tl(n)      r*8  vector of lower bounds on flood threshold
c            tu(n)      r*8  vector of upper bounds on flood threshold
c            alpha      r*8  alpha value to be used in defining plotting
c                              position pp(i) = (i-alpha)/(n+1-2*alpha)
c
c       n.b.  input data are the same as the first 5 arguments to emafit
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       output variables:
c       ------------------------------------------------------------------------
c
c            q(n)       r*8  sorted discharges
c            peX(n)     r*8  estimated exceedance probabilities
c            nt         i*4  number of censoring thresholds
c            thr(nt)    r*8  vector of censoring thresholds
c            peT(nt)    r*8  estimated exceedance probabilities of thresholds
c            nb(nt)     i*4  number of observations censored at each threshold
c
c       n.b.  output data correspond to every observation that is not identified
c             as broken record.  Censored data are assigned relative plotting
c             positions according to their order in the input time series.
c             
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
      subroutine plotposHS(n,ql,qu,tl,tu,alpha,q,peX,nt,thr,peT,nb)

      implicit none
      
      integer
     1 n,i,j,n_true,nmax
     
C JFE      parameter (nmax=10000)
      parameter (nmax=20000)
      
      integer
     1 ix(nmax),nt,nb(nmax)
      
      double precision
     1 ql(n),qu(n),tl(n),tu(n),peX(*),peT(*),infty,
     2 q(nmax),t(nmax),mu,sigma,pp(nmax),xsynth(nmax),
     3 alpha,thr(nmax)
     
      data infty/1.d12/
     
        n_true = 0
      do 10 i=1,n
        if(tl(i) .lt. infty) then
          n_true = n_true + 1
c  for interval data, use average of endpoints to determine order for plotting
          q(n_true) = (ql(i)+qu(i))/2.d0
          t(n_true) = tl(i)
          ix(n_true) = i
        else
          peX(i) = -99.d0
        endif
10    continue

      
      call ppfit2(n_true,q,t,mu,sigma,pp,xsynth,
     1            alpha,nt,thr,peT,nb)
      
      do 20 j=1,n_true
        peX(ix(j)) = 1.d0-pp(j)
20    continue
      
      return
      end
C*_*-*~*_*-*~*             NEW PROGRAM BEGINS HERE            *_*-*~*_*-*~*_*-*~
         SUBROUTINE MDL_MR(N,X,XD,MU,SIGMA,PPALL,XSYNTH,MOMS,QUANTS)

C===============================================================================
C
C        SUBROUTINE TO FIT LN2 DISTRIBUTION TO MULTIPLY-CENSORED DATA
C          USING THE MR METHOD (SEE HELSEL AND COHN, WRR, 1988)
C
C
C        AUTHOR........TIM COHN
C        DATE..........12 FEBRUARY 2003 (TAC)
C
C===============================================================================
C     
C     ***  N.B.  THIS METHOD SHOULD NOT BE APPLIED TO WATER QUALITY DATA 
C     ***          EXCEPT WITH GREAT CARE...  TAC (12 MARCH 2003)
C
C     ***        THE MR/MDL APPROACH ASSUMES A TYPE II CENSORING MODEL;
C     ***          CENSORING RELATED TO WATER QUALITY DATA IS INVARIABLY TYPE I
C
C     ***        SEE WARNINGS OR CONTACT AUTHOR IF THERE IS ANY QUESTION
C     ***          ABOUT THIS.
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        N         I*4       INPUT NUMBER OF OBSERVATIONS (NO-DETECTS COUNT)
C        X(N)      R*8       INPUT VECTOR CONTAINING DATA;  '0' INDICATES 
C                              NO-DETECT
C        XD(N)     R*8       INPUT VECTOR CONTAINING CENSORING THRESHOLD FOR 
C                              EACH PT.
C        MU        R*8       OUTPUT FITTED VALUE OF MU
C        SIGMA     R*8       OUTPUT FITTED VALUE OF SIGMA
C        PPALL(N)  R*8       OUTPUT VECTOR OF PLOTTING POSITIONS CORR. TO X
C                              (SEE HELSEL & COHN, WRR 24(12), 1988)
C        XSYNTH(N) R*8       OUTPUT VECTOR OF REAL AND SYNTHETIC OBSERVATIONS
C                              XSYNTH(I) = X(I); X(I)>=XD(I) 
C                              XSYNTH(I) = EXP(MU+Z(PPALL(I))*SIGMA); X(I)<XD(I)
C        MOMS(2)   R*8       OUTPUT VECTOR OF MEAN(X), VARIANCE(X)
C        QUANTS(5) R*8       OUTPUT VECTOR OF 10TH, 25TH, 50TH, 75TH, AND 90TH
C                              PERCENTILES OF SAMPLE USING WEIBULL PLOTTING POS.
C
C===============================================================================
C
C        N.B.   1)  "SYNTHETIC" VALUES ARE ASSIGNED TO THE "LESS-THAN" VALUES.
C               2)  THE "SYNTHETIC VALUES" MAY EXCEED THE CENSORING THRESHOLD
C                   (THIS OCCURS FOR TWO REASONS:
C                    1)  AN LN2 MODEL IS EMPLOYED TO MODEL THE CENSORED DATA
C                    2)  THE PLOTTING POSITION AND MR FORMULAS ASSUME A
C                        TYPE-II CENSORING MODEL.  IN FACT, CENSORING OF WATER
C                        QUALITY DATA ALMOST INVARIABLY RESULTS FROM A TYPE-I
C                        CENSORING PROCESS)   
C
C===============================================================================

      IMPLICIT NONE

      DOUBLE PRECISION
     1     X(*),XD(*),MU,SIGMA,PPALL(*),XSYNTH(*),MOMS(2),
     2     P(5),QUANTS(5),
     3     XS,XSS
     
      INTEGER
     1     N,I,NP
     
      DATA NP/5/,P/0.10,0.25,0.50,0.75,0.90/
         
      CALL PPFIT(N,X,XD,MU,SIGMA,PPALL,XSYNTH)
         
        XS  = 0.D0
            XSS = 0.D0
          DO 30 I=1,N
            XS  =  XS+XSYNTH(I)
            XSS =  XSS+XSYNTH(I)**2
30    CONTINUE
        MOMS(1)  =  XS/DBLE(N)
        MOMS(2)  =  (XSS-N*MOMS(1)**2)/(N-1.D0)
        
      CALL QTL(XSYNTH,N,P,NP,QUANTS)

      RETURN
      END

         SUBROUTINE PPFIT(N,X,XD,MU,SIGMA,PPALL,XSYNTH)
C===============================================================================
C
C     DUMMY ROUTINE TO REPRODUCE RESULTS FROM PREVIOUS VERSIONS OF PPFIT
C
      
         INTEGER N_T
C JFE
C         PARAMETER (N_T=10000)
         PARAMETER (N_T=20000)

         DOUBLE PRECISION
     1     X(*),XD(*),MU,SIGMA,PPALL(*),XSYNTH(*),PE(N_T),THR(N_T)
     
         INTEGER 
     1     N,NT,NBT(N_T)
         
         CALL PPFIT2(N,X,XD,MU,SIGMA,PPALL,XSYNTH,
     1     0.D0,NT,THR,PE,NBT)
         
         RETURN
         END
C===============================================================================

         SUBROUTINE PPFIT2(N,X,XD,MU,SIGMA,PPALL,XSYNTH,
     1       ALPHA,NT,THR,PE_OFF,NBT)

C===============================================================================
C
C        SUBROUTINE TO FIT LN2 DISTRIBUTION TO MULTIPLY-CENSORED DATA
C
C        AUTHOR........TIM COHN
C        DATE..........APRIL 7, 1986
C        MODIFIED......12 FEBRUARY 2003 (TAC)
C         --REMOVED ERRORS INTRODUCED IN 1986 REVISION TO CODE
C         --IMPLICIT NONE ADDED; ALL VARIABLES DECLARED
C         --DOUBLE PRECISION
C        MODIFIED......10 JUNE 2011 (TAC)
C         --COMPUTE PLOTTING POSITIONS FOR ''GLYPHS''
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        N         I*4       INPUT NUMBER OF OBSERVATIONS (NO-DETECTS COUNT)
C        X(N)      R*8       INPUT VECTOR CONTAINING DATA;  '0' INDICATES
C                              NO-DETECT
C        XD(N)     R*8       INPUT VECTOR CONTAINING CENSORING THRESHOLD FOR 
C                              EACH PT.
C        MU        R*8       OUTPUT FITTED VALUE OF MU
C        SIGMA     R*8       OUTPUT FITTED VALUE OF SIGMA
C        PPALL(N)  R*8       OUTPUT VECTOR OF PLOTTING POSITIONS 
C                              CORRESPONDING TO X
C                              (SEE HELSEL & COHN, WRR 24(12), 1988)
C        XSYNTH(N) R*8       OUTPUT VECTOR OF REAL AND SYNTHETIC OBSERVATIONS
C                              XSYNTH(I) = X(I);  X(I)>=XD(I) 
C                              XSYNTH(I) = EXP(MU+Z(PPALL(I))*SIGMA); X(I)<XD(I)
C        NT        I*4       OUTPUT NUMBER OF DISTINCT THRESHOLD
C        THR(NT)   R*8       OUTPUT VECTOR OF THRESHOLDS
C        PE_OFF(NT)R*8       OUTPUT VECTOR NON-EXCEEDANCE PROBS CORR. TO
C                              THRESHOLDS; NOTE THAT PE(0) HAS BEEN OMITTED
C        NBT(NT)    I*4       OUTPUT NO. NO-DETECTS CENSORED AT THR(I)
C                              (OBSVS SPECIFIED AS "<THR(I)")
C
C===============================================================================
C
C        N.B.   1)  "SYNTHETIC" VALUES ARE ASSIGNED TO THE "LESS-THAN" VALUES.
C               2)  THE "SYNTHETIC VALUES" MAY EXCEED THE CENSORING THRESHOLD
C                   (THIS OCCURS FOR TWO REASONS:
C                    1)  AN LN2 MODEL IS EMPLOYED TO MODEL THE CENSORED DATA
C                    2)  THE PLOTTING POSITION AND MR FORMULAS ASSUME A
C                        TYPE-II CENSORING MODEL.  IN FACT, CENSORING OF WATER
C                        QUALITY DATA ALMOST INVARIABLY RESULTS FROM A TYPE-I
C                        CENSORING PROCESS)   
C
C===============================================================================
C
C        NT        I*4       NUMBER OF THRESHOLDS
C        NBT(NT)   I*4       NUMBER OF OBSERVATIONS BELOW EACH THRESHOLD
C        ND(NT)    I*4       NUMBER OF DETECTS ABOVE THIS THR. AND BELOW NEXT
C        NVAL      I*4       TOTAL NUMBER OF DETECTS
C
C===============================================================================

         IMPLICIT NONE

         INTEGER N_PTS,N_T
C JFE
C         PARAMETER (N_T=10000,N_PTS=10000)
         PARAMETER (N_T=20000,N_PTS=20000)

         DOUBLE PRECISION
     1     X(*),XD(*),MU,SIGMA,
     2     Y(N_PTS),PP(N_PTS),ALPHA,
     3     PPC(N_PTS),PPALL(*),XSYNTH(*),
     4     THR(N_T),PE(0:N_T),PE_OFF(N_T)
     
         INTEGER
     1     N,ND(N_T),NBT(N_T),NB(N_T),NT,NVAL,I

         CALL ARRANGE2(X,XD,N,NT,ND,NB,NVAL,Y,THR,NBT)

         CALL PPLOT2(NT,ND,NB,ALPHA,PP,PPC,PE)
           DO 10 I=1,NT
             PE_OFF(I) = PE(I)      ! ELIMINATE PE(0), WHICH EQUALS 1.0
   10      CONTINUE
         
         CALL PPSOLVE(Y,PP,NVAL,MU,SIGMA)
         
         CALL PLPOS(N,X,XD,PP,PPC,MU,SIGMA,PPALL,XSYNTH)

         RETURN
         END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE ARRANGE2(X,XD,N,NT,ND,NB,NVAL,Y,THRESH,NBT)
C===============================================================================
C
C        SUBROUTINE DESIGNED TO ARRANGE DATA FOR PPFIT, WHICH EMPLOYS DATA TO
C        FIT LN2 DISTRIBUTION TO MULTIPLY-CENSORED DATA
C
C        AUTHOR........TIM COHN
C        DATA..........APRIL 7, 1986
C        MODIFIED......10 JUNE 2011 (TAC)
C         --COMPUTE PLOTTING POSITIONS FOR ''GLYPHS''
C        MODIFIED......12 FEBRUARY 2003 (TAC)
C         --REMOVED ERRORS INTRODUCED IN 1986
C           WHEN CODE WAS REVISED BY DRH.  
C         --IMPLICIT NONE ADDED; ALL VARIABLES DECLARED
C         --DOUBLE PRECISION
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        X(N)      R*8       INPUT VECTOR CONTAINING DATA;  '0' INDICATES 
C                              NO-DETECT
C        XD(N)     R*8       INPUT VECTOR CONTAINING CENSORING THRESHOLD FOR 
C                              EACH PT.
C        N         I*4       INPUT NUMBER OF OBSERVATIONS (NO-DETECTS COUNT)
C        NT        I*4       OUTPUT NUMBER OF THRESHOLDS
C        NB(NT)    I*4       OUTPUT TOTAL NO. OBS. BELOW EACH THRESHOLD
C        NBT(NT)   I*4       OUTPUT NO. OBS. BELOW SPECIFIC THRESHOLD
C        ND(NT)    I*4       OUTPUT NUMBER OF DETECTS ABOVE THIS THR. AND 
C                              BELOW NEXT
C        NVAL      I*4       OUTPUT TOTAL NUMBER OF DETECTS
C        Y(N)      R*8       OUTPUT VECTOR OF SORTED DETECTS
C        THRESH(NT)R*8       OUTPUT VECTOR OF SORTED THRESHOLDS
C
C===============================================================================

         IMPLICIT NONE

         DOUBLE PRECISION
     1      X(*),XD(*),Y(*),
     2      THRESH(200),PLUSINF
     
         INTEGER
C JFE
C     1      N,NT,NVAL,ND(*),NBT(*),NB(*),IP(10000),
     1      N,NT,NVAL,ND(*),NBT(*),NB(*),IP(20000),
     2      I,ICT,J

         DATA PLUSINF/1.0D30/

         CALL PORDER(XD,N,IP)

           ICT          =  0
           J            =  1
           THRESH(1)    =  XD(IP(1))
           NBT(1)       =  0

         DO 10 I=1,N

            IF(XD(IP(I)) .NE. THRESH(J)) THEN
              J    =  J+1
              THRESH(J) =  XD(IP(I))
              ND(J)     =  0
              NBT(J)    =  0
            ENDIF
              
            IF(X(IP(I)) .GE. THRESH(J)) THEN
              ICT       =  ICT+1
              Y(ICT)    =  X(IP(I))
            ELSE
              NBT(J)    =  NBT(J)+1
            ENDIF
   10    CONTINUE

         NT             =  J
         NVAL           =  ICT
         THRESH(NT+1)   =  PLUSINF

         CALL SHELLSORT(Y,ICT)
                
           J       =  1
           ND(1)   =  0
         DO 20 I=1,NVAL
   30      IF(Y(I) .GE. THRESH(J+1)) THEN
              J    =  J+1
              ND(J)=  0
              GOTO 30
           ELSE
              ND(J)=  ND(J)+1
           ENDIF
   20    CONTINUE

           NB(1)       =  NBT(1)
         DO 40 I=2,NT
           NB(I)   =  NB(I-1)+NBT(I)+ND(I-1)
   40    CONTINUE

         RETURN
         END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE PPLOT2(NT,ND,NB,ALPHA,PP,PPC,PE)
C===============================================================================
C
C        SUBROUTINE TO FIND PLOTTING POSITIONS FOR DATA CENSORED AT MULTIPLE 
C          THRESHOLDS
C
C        AUTHOR........TIM COHN
C        DATE..........APRIL 7, 1986
C        MODIFIED......12 FEBRUARY 2003 (TAC)
C         --REMOVED ERRORS INTRODUCED IN 1986
C           WHEN CODE WAS REVISED BY DRH.  
C         --IMPLICIT NONE ADDED; ALL VARIABLES DECLARED
C         --DOUBLE PRECISION
C        MODIFIED......15 JUNE 2011 (TAC)
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        NT        I*4       INPUT NUMBER OF THRESHOLDS
C        ND(NT)    I*4       INPUT NUMBER OF DETECTS ABOVE THIS THR. AND BELOW 
C                              NEXT
C        NB(NT)    I*4       INPUT NUMBER OF OBSERVATIONS BELOW EACH THRESHOLD
C        ALPHA     R*8       INPUT PLOTTING POSITION ALPHA (I-A)/(N+1-2A)
C        PP(NA)    R*8       OUTPUT VECTOR OF PLOTTING POSITIONS
C        PPC(NB)   R*8       OUTPUT VECTOR OF PLOTTING POSITIONS CORRESONDING TO 
C                              X<XD
C        PE(0:NT)  R*8       OUTPUT VECTOR OF PLOTTING POSITIONS CORRESONDING TO 
C                              THRESHOLDS
C
C===============================================================================

         IMPLICIT NONE

         INTEGER N_T
C JFE
C
C         PARAMETER (N_T=10000)
         PARAMETER (N_T=20000)

         DOUBLE PRECISION
     1     PP(*),PE(0:N_T),ALPHA,PD,PPT,PPC(*)
     
         INTEGER
     1     ND(*),NB(*),NT,
     2     ISUM,I,J,JSUM,NBT


           PE(NT+1)  =  0.D0
         DO 10 I=NT,1,-1
           PD        =  ND(I)/FLOAT(MAX(1,ND(I))+NB(I))
           PE(I)     =  PE(I+1) + (1.D0-PE(I+1)) * PD
   10    CONTINUE

           ISUM    =  0
           JSUM    =  0
           PE(0)   =  1.00D0

         DO 20 I=1,NT
           DO 30 J=1,ND(I)
               PPT  =  (J-ALPHA)/(ND(I) + 1.0 - 2.0*ALPHA)
               PP(ISUM+J)    =  (1.D0-PE(I)) + (PE(I)-PE(I+1))*PPT
   30      CONTINUE
               ISUM =  ISUM+ND(I)

           IF(I .EQ. 1) THEN
              NBT = NB(1)
           ELSE
              NBT = NB(I)-NB(I-1)-ND(I-1)
           ENDIF
           
           DO 40 J=1,NBT
             PPT  = (J-ALPHA)/(NBT+1.0-2*ALPHA)
             PPC(JSUM+J) = (1-PE(I))*PPT
   40      CONTINUE
             JSUM=JSUM+NBT

   20    CONTINUE

         RETURN
         END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE PPSOLVE(Y,PP,NVAL,MU,SIGMA)
C===============================================================================
C
C        SUBROUTINE TO SOLVE FOR LEAST SQUARES ESTIMATES OF MU AND SIGMA
C        ASSUMING A LINEAR MODEL
C
C        AUTHOR........TIM COHN
C        DATE..........APRIL 7, 1986
C        MODIFIED......12 FEBRUARY 2003 (TAC)
C         --REMOVED ERRORS INTRODUCED IN 1986
C           WHEN CODE WAS REVISED BY DRH.  
C         --IMPLICIT NONE ADDED; ALL VARIABLES DECLARED
C         --DOUBLE PRECISION
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        Y(N)      R*8       INPUT VECTOR CONTAINING SORTED DETECT-DATA
C        PP(NA)    R*8       INPUT VECTOR OF PLOTTING POSITIONS
C        NVAL      I*4       INPUT TOTAL NUMBER OF DETECTS
C        MU        R*8       OUTPUT FITTED VALUE OF MU
C        SIGMA     R*8       OUTPUT FITTED VALUE OF SIGMA
C
C===============================================================================

         IMPLICIT NONE

         INTEGER N_PTS
C JFE
C         PARAMETER (N_PTS=10000)
         PARAMETER (N_PTS=20000)

         DOUBLE PRECISION
     1     Y(*),PP(*),MU,SIGMA,
     2     PSIG(N_PTS),LOGY(N_PTS),
     3     FP_Z_ICDF,
     4     SUMY,SUMP,SUMPY,SUMP2,YBAR,PBAR
     
         INTEGER
     1     I,NVAL

              SUMY =  0.D0
              SUMP =  0.D0
              SUMPY=  0.D0
              SUMP2=  0.D0

         DO 10 I=1,NVAL
              PSIG(I)   =  FP_Z_ICDF(PP(I))
              LOGY(I)   =  LOG(Y(I))
              SUMY      =  SUMY    +  LOGY(I)
              SUMPY     =  SUMPY   +  LOGY(I)*PSIG(I)
              SUMP      =  SUMP    +  PSIG(I)
              SUMP2     =  SUMP2   +  PSIG(I)**2
   10    CONTINUE
              YBAR      =  SUMY/NVAL
              PBAR      =  SUMP/NVAL

              SIGMA     =  (SUMPY-NVAL*PBAR*YBAR)/(SUMP2-NVAL*PBAR**2)
              MU        =  YBAR - SIGMA*PBAR
         RETURN
         END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE PLPOS(N,X,XD,PP,PPC,MU,SIGMA,PPALL,XSYNTH)
C===============================================================================
C
C        SUBROUTINE PLPOS GIVES THE PLOTTING POSITIONS FOR ALL OF THE
C          OBSERVATIONS, CENSORED AND UNCENSORED, IN THE ORIGINAL INPUT
C          DATA SET
C
C
C        AUTHOR....TIM COHN
C        DATE......12 FEBRUARY 2003 (TAC)
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        N         I*4       INPUT NUMBER OF OBSERVATIONS (NO-DETECTS COUNT)
C        X(N)      R*8       INPUT VECTOR CONTAINING DATA;  '0' INDICATES 
C                              NO-DETECT
C        XD(N)     R*8       INPUT VECTOR CONTAINING CENSORING THRESHOLD FOR 
C                              EACH PT.
C        PP(NA)    R*8       INPUT VECTOR OF PLOTTING POSITIONS CORRESONDING TO 
C                              X>=XD
C                              (SEE HELSEL & COHN, WRR 24(12), 1988)
C        PPC(NB)   R*8       INPUT VECTOR OF PLOTTING POSITIONS CORRESONDING TO 
C                              X<XD
C        MU        R*8       INPUT FITTED VALUE OF MU
C        SIGMA     R*8       INPUT FITTED VALUE OF SIGMA
C        PPALL(N)  R*8       OUTPUT VECTOR OF PLOTTING POSITIONS CORRESONDING TO X
C                              (SEE HELSEL & COHN, WRR 24(12), 1988)
C        XSYNTH    R*8       OUTPUT VECTOR OF REAL AND SYNTHETIC OBSERVATIONS
C                              XSYNTH(I) = X(I);                       
C                              X(I)>=XD(I) 
C                              XSYNTH(I) = EXP(MU+Z(PPALL(I))*SIGMA);  
C                              X(I)<XD(I)
C
C===============================================================================

         IMPLICIT NONE
         
         INTEGER N_PTS
C JFE
C
C         PARAMETER (N_PTS=10000)
         PARAMETER (N_PTS=20000)

         DOUBLE PRECISION 
     1     X(*),XD(*),MU,SIGMA,PPALL(*),PP(*),PPC(*),XSYNTH(*),
     2     FP_Z_ICDF
         
         INTEGER
     1     IX(N_PTS),IXD(N_PTS),N,NA,NB,I

         CALL PORDER(X,N,IX)
         CALL PORDER(XD,N,IXD)
         
           NA  =  0
           NB  =  0
         DO 10 I=1,N
           IF(X(IX(I)) .GE. XD(IX(I))) THEN
             NA            =  NA+1
             PPALL(IX(I))  =  PP(NA)
             XSYNTH(IX(I)) =  X(IX(I))
           ENDIF
           IF(X(IXD(I)) .LT. XD(IXD(I))) THEN
             NB            =  NB+1
             PPALL(IXD(I)) =  PPC(NB)
             XSYNTH(IXD(I))=  EXP(MU+FP_Z_ICDF(PPC(NB))*SIGMA)
           ENDIF
10       CONTINUE

         RETURN
         END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE PORDER(X,N,IX)
C===============================================================================
C
C        SUBROUTINE PORDER GIVES THE PERMUTATION OF A USER-SUPPLIED
C        VECTOR.  TREE-SORT ALGORITHM IS EMPLOYED (WHY NOT?)
C        RETURN RESULTING IX ARRAY WITH THE ASCENDING PERMUTATIONS
C
C          X(IX(1)) = SMALLEST VALUE IN X
C          X(IX(N)) = LARGEST VALUE IN X
C
C
C        AUTHOR....TIM COHN
C        DATE......APRIL 1,  1986
C        REVISED...AUGUST 9, 1986      (TAC)
C        MODIFIED......12 FEBRUARY 2003 (TAC)
C         --REMOVED ERRORS INTRODUCED IN 1986
C           WHEN CODE WAS REVISED BY DRH.  
C         --IMPLICIT NONE ADDED; ALL VARIABLES DECLARED
C         --DOUBLE PRECISION
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        X(N)      R*8       INPUT VECTOR OF UNORDERED DATA
C        N         I*4       INPUT NUMBER OF OBSERVATIONS IN X
C        IX(N)     I*4       OUTPUT VECTOR CONTAINING PERMUTATION OF X
C                               X(IX(1)) IS THE SMALLEST VALUE OF X
C                                 . . . 
C                               X(IX(N)) IS THE LARGEST VALUE OF X
C
C===============================================================================

         IMPLICIT NONE
         
         INTEGER N_PTS
C JFE
C
C         PARAMETER (N_PTS=10000)
         PARAMETER (N_PTS=20000)

         DOUBLE PRECISION X(*)
         
         INTEGER
     1     L(0:N_PTS),R(0:N_PTS),P(0:N_PTS),
     2     IX(*),I2,I,INDX,N,ICT

C===============================================================================
C
C    FIRST CHECK TO SEE IF WE HAVE AN ORDERED DATA VECTOR TO BEGIN WITH
C
         DO 50 I2=2,N
              IX(I2)    =  I2
              IF(X(I2) .LT. X(I2-1)) GOTO 1
   50    CONTINUE
              IX(1)     =  1
              RETURN

    1    CONTINUE

         L(1) =  0
         R(1) =  0
         P(1) =  0

        DO 10 I=2,N
         INDX =  1
         L(I)   =  0
         R(I)   =  0

   20    CONTINUE
         IF(X(I) .GE. X(INDX)) THEN
              IF(R(INDX) .EQ. 0) THEN
                   R(INDX)   =  I
                   P(I)      =  INDX
                   GOTO 10
              ELSE
                   INDX =  R(INDX)
                   GOTO 20
              ENDIF
         ELSE
              IF(L(INDX) .EQ. 0) THEN
                   L(INDX)   =  I
                   P(I)      =  INDX
                   GOTO 10
              ELSE
                   INDX =  L(INDX)
                   GOTO 20
              ENDIF
         ENDIF
   10   CONTINUE

          INDX =  1
         DO 40 ICT=1,N

   30    CONTINUE
         IF(L(INDX) .EQ. 0) THEN
              IX(ICT)     =  INDX
              P(R(INDX))   =  P(INDX)
              L(P(INDX))   =  R(INDX)
              INDX         =  P(INDX)
         ELSE
              INDX =  L(INDX)
              GOTO 30
         ENDIF
   40    CONTINUE
         RETURN
         END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE SHELLSORT(X,N)
C===============================================================================
C
C        SUBROUTINE TO SORT A VECTOR IN PLACE
C
C        AUTHOR........TIM COHN
C        DATE..........12 FEBRUARY 2003 (TAC)
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        X(N)      R*8       INPUT/OUTPUT VECTOR OF DATA
C        N         I*4       INPUT SAMPLE SIZE
C
C===============================================================================

      IMPLICIT NONE
      
      DOUBLE PRECISION X(*),R
      
      INTEGER N,M,I,J,S,Z
      
      M=N
  100 M=INT(M/2)
      IF (M.EQ.0) GOTO 210
      DO 190 S=1,M
  130 I=S
      J=S+M
      Z=0
  140 IF (X(I).LE.X(J)) GOTO 160
      Z=1
      R=X(I)
      X(I)=X(J)
      X(J)=R
  160 I=J
      J=J+M
      IF (J.LT.(N+1)) GOTO 140
      IF (Z.EQ.1) GOTO 130
  190 CONTINUE
      GOTO 100
  210 RETURN
      END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
         SUBROUTINE QTL(X,N,P,NP,Q)
C===============================================================================
C
C        SUBROUTINE TO COMPUTE QUANTILES OF A VECTOR
C        THIS EMPLOYS THE WEIBULL OR GUMBEL PLOTTING POSITIONS (I/(N+1))
C
C        AUTHOR........TIM COHN
C        DATE..........12 FEBRUARY 2003 (TAC)
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)   
C
C
C===============================================================================
C
C        X(N)      R*8       INPUT VECTOR OF DATA
C        N         I*4       INPUT SAMPLE SIZE
C        P(N)      R*8       INPUT VECTOR OF PROBABILITIES (PERCENTILES/100) TO 
C                              ESTIMATE
C        Q(N)      R*8       OUTPUT VECTOR OF QUANTILES 
C
C===============================================================================
C
C        N.B.   1)  WEIBULL (AKA "GUMBEL") PLOTTING POSTIONS ARE USED, WHICH 
C                   CORRESPOND TO ALPHA=0.0 IN THE GENERAL EQUATION FOR 
C                   PLOTTING POSTIONS PP(I) =(I-ALPHA)/(N+1-2*ALPHA)
C                   EXCEL, AND MANY OTHER STATISTICAL PACKAGES, EMPLOY THE 
C                   "HAZEN" PLOTTING POSITIONS (ALPHA=0.5).
C                   
C===============================================================================

         IMPLICIT NONE
         
         INTEGER N_PTS
C JFE
C
C         PARAMETER (N_PTS=10000)
         PARAMETER (N_PTS=20000)

         DOUBLE PRECISION 
     1     X(*),P(*),Q(*),XS(0:N_PTS),FN,ALPHA
     
         INTEGER 
     1     I,INDX,NP,N
     
         CALL DSVRGN(N,X,XS(1))
           XS(0)     = -9.0D99
           XS(N+1)   =  9.0D99
         
         DO 20 I=1,NP
           FN        =  MAX(0.D0,MIN(1.D0,P(I)))*(N+1.D0)
           INDX      =  INT(FN)
           ALPHA     =  FN-INDX
           Q(I)      =  (1.D0-ALPHA)*XS(INDX)+ALPHA*XS(INDX+1)
20       CONTINUE

         RETURN
         END
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    silly little function to print asterisks for p-value significance
c
c    timothy a. cohn        9 mar 2011
c
      character*4 function signif(p)
      double precision p
      if(p .lt. 0.01) then
        signif = ' ***'
      else if(p .lt. 0.05d0) then
        signif = ' ** '
      else if(p .lt. 0.10d0) then
        signif = ' *  '
      else if(p .ge. 0.10d0) then
        signif = '    '
      endif
        return
      end
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    another silly little function to grab last 4 characters in filename
c      (usually the ".ext") and make it available.
c
c    timothy a. cohn        9 mar 2011
c
      character*4 function filext(fname)
      character*(*) fname
      integer ichar
        ichar=len_trim(fname)
      filext=fname(max(1,(ichar-3)):ichar)
      return
      end
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c****|subroutine EDFinit
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c     Program to estimate probabilities and quantiles of distributions
c       based on finite samples
c    
c       
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    development history
c
c    timothy a. cohn        15 jul 2011 
c       modified            18 jul 2011
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       input variables:
c       ------------------------------------------------------------------------
c            n          i*4  number of observations (censored, uncensored, or 
c                              other)
c            x(n)       r*8  vector of edf data
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       output variables:
c       ------------------------------------------------------------------------
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
      subroutine EDFinit(n,x)

      implicit none

      integer
     1  n,                                                ! input variables
     2  pn,i,no
     
      logical initc

      parameter (pn=1000000)
      
      double precision 
     1  x(*),xc,pc,pnep
     
      save
     
      common /EDF001/no,initc,xc(0:pn),pc(0:pn),pnep(0:pn)
      
c  initiate arrays for work

      if(n .lt. 1) then
        write(*,*) 'EDFinit: Need more than 0 observations'
        stop
      else if(n .lt. 1) then
        write(*,*) 'EDFinit: Maximum Number of Datapoints: ',n
        stop
      endif
      
      no    = n
      call dsvrgn(no,x,xc(1))
      
      do 10 i=1,no
          pc(i)  = dble(i)/dble(no) ! pc
10    continue

      do 20 i=no-1,1,-1
        if( xc(i) .lt. xc(i+1) ) then
          pnep(i) = pc(i)        ! pnep(i) = non-exceedance prob of xc(i)
        else
          pnep(i) = pnep(i+1)   ! dealing with ties
        endif
20    continue
        pnep(no) = 1.d0
      
      call srand(123457)      ! initialize random number generator with new seed
      initc = .TRUE.
      
      return
      end
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
      blockdata EDFblk
      implicit none

      integer pn,no
     
      parameter (pn=1000000)
      
      double precision xc,pc,pnep
     
      logical initc
     
      common /EDF001/no,initc,xc(0:pn),pc(0:pn),pnep(0:pn)
      
      data initc/.FALSE./,xc(0)/-1.d99/,pc(0)/0.d0/,pnep(0)/0.d0/
      
      end
      
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c****|double precision function pEDF
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c     Program to estimate probabilities of distributions
c       based on finite samples
c    
c       
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    development history
c
c    timothy a. cohn        15 jul 2011 
c       modified            18 jul 2011
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       input variables:
c       ------------------------------------------------------------------------
c            x0         r*8  x value for which F(x0) is to be computed
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       output
c       ------------------------------------------------------------------------
c            pEDF    r*8  F(x0)
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
      double precision function pEDF(x0)
      
      implicit none

      logical initc
     
      integer
     1  pn,i,no,itop,ibot,iguess
     
      parameter (pn=1000000)
      
      double precision 
     1  x0,xc,result,pc,pnep
     
      save
     
      common /EDF001/no,initc,xc(0:pn),pc(0:pn),pnep(0:pn)
      
      if(.not. initc) then
        write(*,*) 'pEDF: EDF needs to be initiated before use'
        stop
      endif
      
        ibot = 0
        itop = no
      do 10 i=1,1000
        if(itop-ibot .le. 1) then
          result = pc(itop)
          goto 99
        endif
          iguess = (ibot+itop)/2
        if(x0 .gt. xc(iguess)) then
          ibot = iguess
        else if(x0 .lt. xc(iguess)) then
          itop = iguess
        else
          result = pnep(iguess)
          goto 99
        endif
10        continue
99    continue
        pEDF = result
      return
      end
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c****|double precision function qEDF
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c     Program to estimate quantiles of distributions
c       based on finite samples
c    
c       
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    development history
c
c    timothy a. cohn        15 jul 2011 
c       modified            18 jul 2011
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       input variables:
c       ------------------------------------------------------------------------
c            p0         r*8  x value for which F(x0) is to be computed
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       output
c       ------------------------------------------------------------------------
c            qEDF   r*8  F(x0)
c
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c   N.B. Quantile function corresponds to R quantile function
c                  quantile(xc,p0,type=4)
c        The quantiles are interpolated between datapoints 
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
      double precision function qEDF(p0)
      
      implicit none

      integer
     1  pn,i,no,itop,ibot,iguess
     
      logical initc
     
      parameter (pn=1000000)
      
      double precision 
     1  p0,xc,pc,pnep,result
     
      save
     
      common /EDF001/no,initc,xc(0:pn),pc(0:pn),pnep(0:pn)
      
      if(.not. initc) then
        write(*,*) 'qEDF: EDF needs to be initiated before use'
        stop
      endif
      
        ibot = 0
        itop = no
      do 10 i=1,1000
        if(itop-ibot .le. 1) then
          if(pc(itop) .ne. pc(ibot)) then
            result = ((p0-pc(ibot))*xc(itop)+(pc(itop)-p0)*xc(ibot))/
     1                          (pc(itop)-pc(ibot))
          else
            result = pc(ibot)
          endif
          goto 99
        endif
          iguess = (ibot+itop)/2
        if(p0 .gt. pc(iguess)) then
          ibot = iguess
        else if(p0 .lt. pc(iguess)) then
          itop = iguess
        else
          result = xc(iguess)
          goto 99
        endif
10        continue
99    continue
        qEDF = result
      return
      end
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c****|double precision function rEDF
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c     Program to generate random samples (sampled with replacement)
c       from EDF  
c       
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c    development history
c
c    timothy a. cohn        15 jul 2011 
c       modified            18 jul 2011
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
c       input variables:
c       ------------------------------------------------------------------------
c
c
c****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
c
      double precision function rEDF()
      
      implicit none

      integer
     1  no,pn,indx
     
      logical initc
     
      parameter (pn=1000000)
      
      double precision 
     1  xc,pc,pnep
     
      real*4 rand
     
      save
     
      common /EDF001/no,initc,xc(0:pn),pc(0:pn),pnep(0:pn)

      if(.not. initc) then
        write(*,*) 'rEDF: EDF needs to be initiated before use'
        stop
      endif

        indx = no*rand() + 1
        rEDF = xc(indx)
      return
      end
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
       SUBROUTINE NORMQUAD(N,T,W)
C===============================================================================
C
C        SUBROUTINE NORMQUAD COMPUTES QUADRATURE NODES AND WEIGHTS 
C          CORRESPONDING TO THE STANDARD NORMAL DISTRIBUTION
C
C          INT(-INF,INF) g(z)f(z)dz = sum(i=1,n) w(i)*g(t(i))
C
C
C        AUTHOR....TIM COHN
C        DATE......05 OCTOBER 2012
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        N         I*4       INPUT NUMBER OF NODES/ORDER OF OTHOGONAL POLYNOMIAL
C        T(N)      R*8       OUTPUT VECTOR OF NODES
C        W(N)      R*8       OUTPUT VECTOR OF WEIGHTS
C
C===============================================================================
C           THIS SET OF ROUTINES COMPUTES THE NODES T(J) AND WEIGHTS
C        W(J) FOR GAUSSIAN-TYPE QUADRATURE RULES WITH PRE-ASSIGNED
C        NODES.  THESE ARE USED WHEN ONE WISHES TO APPROXIMATE
C
C                 INTEGRAL (FROM A TO B)  F(X) W(X) DX
C
C                              N
C        BY                   SUM W  F(T )
C                             J=1  J    J
C
C        (NOTE W(X) AND W(J) HAVE NO CONNECTION WITH EACH OTHER.)
C        HERE W(X) IS ONE OF SIX POSSIBLE NON-NEGATIVE WEIGHT
C        FUNCTIONS (LISTED BELOW), AND F(X) IS THE
C        FUNCTION TO BE INTEGRATED.  GAUSSIAN QUADRATURE IS PARTICULARLY
C        USEFUL ON INFINITE INTERVALS (WITH APPROPRIATE WEIGHT
C        FUNCTIONS), SINCE THEN OTHER TECHNIQUES OFTEN FAIL.
C
C           ASSOCIATED WITH EACH WEIGHT FUNCTION W(X) IS A SET OF
C        ORTHOGONAL POLYNOMIALS.  THE NODES T(J) ARE JUST THE ZEROES
C        OF THE PROPER N-TH DEGREE POLYNOMIAL.
C
C     INPUT PARAMETERS (ALL REAL NUMBERS ARE IN DOUBLE PRECISION)
C
C        KIND     AN INTEGER BETWEEN 1 AND 6 GIVING THE TYPE OF
C                 QUADRATURE RULE:
C
C        KIND = 1:  LEGENDRE QUADRATURE, W(X) = 1 ON (-1, 1)
C        KIND = 2:  CHEBYSHEV QUADRATURE OF THE FIRST KIND
C                   W(X) = 1/SQRT(1 - X*X) ON (-1, +1)
C        KIND = 3:  CHEBYSHEV QUADRATURE OF THE SECOND KIND
C                   W(X) = SQRT(1 - X*X) ON (-1, 1)
C        KIND = 4:  HERMITE QUADRATURE, W(X) = EXP(-X*X) ON
C                   (-INFINITY, +INFINITY)
C        KIND = 5:  JACOBI QUADRATURE, W(X) = (1-X)**ALPHA * (1+X)**
C                   BETA ON (-1, 1), ALPHA, BETA .GT. -1.
C                   NOTE: KIND=2 AND 3 ARE A SPECIAL CASE OF THIS.
C        KIND = 6:  GENERALIZED LAGUERRE QUADRATURE, W(X) = EXP(-X)*
C                   X**ALPHA ON (0, +INFINITY), ALPHA .GT. -1
C===============================================================================
C       
       IMPLICIT NONE
       
       DOUBLE PRECISION T(*),W(*),ENDPTS(2),B(1000),SQ2,SQPI
       INTEGER I,KIND,KPTS,N
       
       DATA KPTS/0/,ENDPTS/2*0.D0/
       DATA SQ2/1.41421356237/,SQPI/0.56418958354/
       
         KIND=4
       CALL GAUSSQ(KIND, N, 0.D0, 0.D0, KPTS, ENDPTS, B, T, W)
       DO 10 I=1,N
         T(I) = T(I)*SQ2
         W(I) = W(I)*SQPI
10     CONTINUE
       
       RETURN
       END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
       SUBROUTINE GAMMAQUAD(N,ALPHA,BETA,T,W)
C===============================================================================
C
C        SUBROUTINE GAMMAQUAD COMPUTES QUADRATURE NODES AND WEIGHTS 
C          CORRESPONDING TO A GAMMA DISTRIBUTION WITH PARAMETERS
C          ALPHA, BETA
C
C          INT(0,INF) g(y)f(y)dy = sum(i=1,n) w(i)*g(t(i))
C
C
C        AUTHOR....TIM COHN
C        DATE......05 OCTOBER 2012
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        N         I*4       INPUT NUMBER OF NODES/ORDER OF OTHOGONAL POLYNOMIAL
C        ALPHA     R*8       SHAPE PARAMETER OF GAMMA DISTRIBUTION
C        BETA      R*8       SCALE PARAMETER OF GAMMA DISTRIBUTION
C        T(N)      R*8       OUTPUT VECTOR OF NODES
C        W(N)      R*8       OUTPUT VECTOR OF WEIGHTS
C
C===============================================================================
C       
       IMPLICIT NONE
       
       DOUBLE PRECISION 
     1   T(*),W(*),ALPHA,BETA,ENDPTS(2),B2(1000),SQ2,SQPI,DGI,DGAMMA,
     2   ALPHAMAX,A,B,C
       INTEGER I,KIND,KPTS,N
       DATA ALPHAMAX/160.D0/
       
       DATA KPTS/0/,ENDPTS/2*0.D0/
       DATA SQ2/1.41421356237/,SQPI/0.56418958354/
       
       IF(ALPHA .LE. ALPHAMAX) THEN
         A = ALPHA
         B = BETA
         C = 0.D0
       ELSE
         A = ALPHAMAX
         B = BETA*SQRT(ALPHA/A)
         C = ALPHA*BETA-A*B
       ENDIF
           KIND=6
         CALL GAUSSQ(KIND, N, A-1.D0, 0.D0, KPTS, ENDPTS, B2, T, W)
           DGI = 1.D0/DGAMMA(A)
         DO 10 I=1,N
           T(I) = C + T(I)*B
           W(I) = W(I)*DGI
10       CONTINUE
       
       RETURN
       END
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
       SUBROUTINE CHOL33(S,V,IFLAG)
C===============================================================================
C
C        SUBROUTINE CHOL33 COMPUTES A SORT-OF CHOLESKY DECOMPOSITION (NOT LU)
C          OF A SYMMETRIC POSITIVE-DEFINITE 3x3 MATRIX WITH THE FOLLOWING
C          PROPERTIES:
C            
C          ASSUME WE START WITH A (Nx3) VARIABLE X (THINK {M,S2,G})
C          X HAS 3x3 COVARIANCE MATRIX S AND P=CHOL33(S,P)
C            --THE COLUMNS OF XxP ARE ORTHOGONAL
c            --THE SECOND COLUMN OF P HAS TWO ZEROES (CHOLESKY ITS COL 1)
C
C        N.B. THIS ROUTINE IS USEFUL FOR ONLY ONE PURPOSE: GENERATING
C             SAMPLES USING QUADRATURE POINTS WHERE THE MIDDLE VARIABLE
C             IS GAMMA AND THE TWO OUTER ONES ARE BOTH NORMAL
C             NO CHECK IS DONE TO VERIFY THAT MATRIX HAS NON-ZERO
C             MAIN DIAGONAL OR ANYTHING ELSE...
C
C        FOR THE ALGEBRA, SEE COHN [2012] (NOTEBOOK 3)
C          
C        AUTHOR....TIM COHN
C        DATE......05 OCTOBER 2012
C
C===============================================================================
C
C     PROPERTY OF US GOVERNMENT, U.S. GEOLOGICAL SURVEY
C
C     *** DO NOT MODIFY WITHOUT AUTHOR''S CONSENT ***
C
C     AUTHOR CAN BE CONTACTED AT:  TACOHN@USGS.GOV (703/648-5711)     
C
C===============================================================================
C
C        S(3,3)    R*8       INPUT MATRIX
C        V(3,3)    R*8       OUTPUT MATRIX, MODIFIED CHOLESKY DECOMPOSITION OF S
C        IFLAG     I*4       OUTPUT FLAG; 0 IF OK; 1 IF S NOT PSD
C
C  THE OUTPUT IS NOT UPPER TRIANGULAR, BUT RATHER UT ONCE COLUMNS 1 AND 2 ARE 
C  SWAPPED AND ROWS 1 AND 2 ARE SWAPPED. 
C
C		> chol(var(x))
C				 [,1]      [,2]        [,3]
C		[1,] 0.938117 0.1361061 -0.04404994
C		[2,] 0.000000 0.9575892 -0.07947178
C		[3,] 0.000000 0.0000000  0.84128956
C		> chol33(var(x))
C				  [,1]      [,2]        [,3]
C		[1,] 0.9287822 0.0000000 -0.03242836
C		[2,] 0.1320116 0.9672135 -0.08487969
C		[3,] 0.0000000 0.0000000  0.84128956
C
C===============================================================================
C
      IMPLICIT NONE
      DOUBLE PRECISION S(3,3),V(3,3),T1
      INTEGER IFLAG
      
      IF(S(1,1).LE.0.D0 .OR. S(2,2).LE.0.D0 .OR. S(3,3).LE.0.D0) GOTO 99
      
      V(1,2)=0.D0
      V(3,2)=0.D0
      V(3,1)=0.D0
      V(2,2)=SQRT(S(2,2))
      V(2,1)=S(2,1)/V(2,2)
      V(2,3)=S(2,3)/V(2,2)
        T1=S(1,1)-V(2,1)**2
        IF(T1 .LT. 0.D0) GOTO 99
      V(1,1)=SQRT(T1)
      V(1,3)=(S(3,1)-V(2,3)*V(2,1))/V(1,1)
        T1=S(3,3)-V(2,3)**2-V(1,3)**2
        IF(T1 .LT. 0.D0) GOTO 99
      V(3,3)=SQRT(T1)
        IFLAG=0
      RETURN
99    CONTINUE
C       WRITE(*,*) 'MATRIX S NOT POS-SEMI-D; (CHOL33)'
        IFLAG=1 
        RETURN
      END
C****
      DOUBLE PRECISION FUNCTION MEANW (N,X,W)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C  COMPUTES THE WEIGHTED MEAN OF A VECTOR X
C
      IMPLICIT NONE
	  INTEGER N,I
	  DOUBLE PRECISION X(*),W(*),WSUM,XSUM

        WSUM = 0.D0
        XSUM = 0.D0
      DO 10 I=1,N
        WSUM = WSUM+W(I)
        XSUM = XSUM+W(I)*X(I)
10    CONTINUE
        MEANW=XSUM/WSUM
      RETURN
	  END
C****
      DOUBLE PRECISION FUNCTION COVW (N,X,Y,W)
C****|===|====-====|====-====|====-====|====-====|====-====|====-====|==////////
C
C  COMPUTES THE WEIGHTED MEAN OF A VECTOR X
C
      IMPLICIT NONE
	  INTEGER N,I
	  DOUBLE PRECISION X(*),Y(*),W(*),WSUM,XSUM,MEANW,MX,MY

        MX  = MEANW(N,X,W)
        MY  = MEANW(N,Y,W)      
        WSUM = 0.D0
        XSUM = 0.D0
      DO 10 I=1,N
        WSUM = WSUM+W(I)
        XSUM = XSUM+W(I)*(X(I)-MX)*(Y(I)-MY)
10    CONTINUE
        COVW = XSUM/WSUM
      RETURN
	  END