SUBROUTINE VDTINI (
*
*  inputs
*
     :                          X, Y, NPTS,
*
*  outputs
*
     :                          COEFF, NCOEFF, NVAR, STATUS)
*
*  Module Number: 15.5.2.2
*
*  Module Name: deadtimev
*
*  Keyphrase:
*  ----------
*  obtain an initial estimate of dead time for the least square operation
*
*  Description:
*  ------------
*  Estimate the initial value of dead time and the ratio between digital
*  count rate and analog count rate.
*  The pair pulse correction formula used here is:
*
*       y = ratio * x / (1 + ratio * x * tau)
*
*  where x is the "true" count rate, y is the observed countrate, tau is the
*  dead time (=coeff(1)), and ratio is the ratio between the digital and analog
*  count rates at low count rate (=coeff(2)).
*
*  FORTRAN Name: VDTINI.FOR
*
*  Keywords of Accessed Files and Tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       None
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       None
*  SDAS:
*       UMSPUT
*  OTHERS:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       11-15-85     J.-C. Hsu        Design and coding
*     2       07-20-86     J.-C. Hsu        add ratio as one parameter
*     3       08-20-87     J.-C. Hsu        F77 standard
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --number of input data points
        INTEGER                 NPTS
*                                       --analog count rate array
        DOUBLE PRECISION        X(1),
*                                       --digital count rate array
     :                          Y(1)
*
*== output:
*                                       --number of coefficients
        INTEGER                 NCOEFF,
*                                       --number of variables
     :                          NVAR,
*                                       --error status
     :                          STATUS
*                                       --coefficients
        DOUBLE PRECISION        COEFF(*)
*
*== local:
*                                       --number of high count rate points
        INTEGER                 NHIGH,
     :                          INDX(2000),
     :                          STATOK,
*                                       --loop indices
     :                          I, J
*                                       --count rate above which is regarded as
*                                       --high count rate
        DOUBLE PRECISION        HICT,
*                                       --minimum of x and y arrays
     :                          XMIN, YMIN,
     :                          SUM
*                                       --error message context
        CHARACTER*130           CONTXT, MESS
*=========================begin hsp.inc=========================================
*                                       --status return code
        INTEGER         OK, ERRNUM(20)
*                                       --message destination and priority
        INTEGER         DEST, PRIO
        DATA OK /0/
        DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
     :               711, 712, 713, 714, 715, 716, 717, 718, 719, 720/
        DATA DEST, PRIO /1, 0/
*=========================end hsp.inc===========================================
*-------------------------------------------------------------------------------
*
*  initialize the parameters
*
        NCOEFF = 2
        NVAR = 2
        HICT = 2.D5
        NHIGH = 0
        SUM = 0.
*
*  check if there is enough data points
*
        IF (NPTS .LE. NCOEFF) THEN
            STATUS = ERRNUM(1)
            CONTXT = 'not enough data points to perform least square '
     :                  // 'fitting'
            GO TO 999
        END IF
*
*  estimate the ratio from the lowest count rate
*
        XMIN = X(1)
        DO 10 I = 1, NPTS
            IF (X(I) .LT. XMIN) THEN
                XMIN = X(I)
                YMIN = Y(I)
            END IF
   10   CONTINUE
        COEFF(2) = YMIN / XMIN
*
*  estimate dead time from high count rate data
*
        DO 20 I = 1, NPTS
            IF (X(I) .GT. HICT) THEN
                NHIGH = NHIGH + 1
                INDX(NHIGH) = I
                SUM = SUM + (COEFF(2) * X(I) - Y(I)) /
     :                  (COEFF(2) * X(I) * Y(I))
            END IF
   20   CONTINUE
*
*  check if there is any high count rate data
*
        IF (NHIGH .EQ. 0) THEN
            STATUS = ERRNUM(2)
            CONTXT = 'observed data all have too low countrates, '//
     :                  'unsuitable for dead time calibration'
            GO TO 999
        END IF
*
*  iterate one more time to get a better initial guess
*
        COEFF(1) = SUM / DBLE (NHIGH)
        COEFF(2) = YMIN / (XMIN * (1. - YMIN * COEFF(1)))
*
        SUM = 0.
*
        NHIGH = MIN (NHIGH, 2000)
        DO 30 I = 1, NHIGH
            J = INDX(I)
            SUM = SUM + (COEFF(2) * X(J) - Y(J)) /
     :                  (COEFF(2) * X(J) * Y(J))
   30   CONTINUE
        COEFF(1) = SUM / DBLE(NHIGH)
*
        STATUS = OK
        GO TO 1000
*
*  write error message
*
  999   MESS = 'VDTINI: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,     STATOK)
*
 1000   RETURN
        END