SUBROUTINE VPHAIN (
*
*  inputs
*
     :                          X, Y, NPTS,
*
*  outputs
*
     :                          COEFF, NCOEFF, NVAR, STATUS)
*
*  Module Number: 15.3.1.2
*
*  Module Name: phav
*
*  Keyphrase:
*  ----------
*  To obtain initial values of the coefficients for the least square
*  procedure to obtain the pulse height distribution function
*
*  Description:
*  ------------
*  Figure out the initial values of the coefficients which specify the
*  integrated pulse height distribution (IPHD).  The pulse height distribution
*  (PHD) function is assumed to be:
*
*       PHD = P/R * exp-((x-Q)**2/2*R**2) + S/T * exp(-x/T)
*
*  where P = coeff(1), Q = coeff(2)..., T = coeff(5).
*  Furthermore, a constant term coef(6) is added to the INTEGRATED PHD.
*  Some of the initial values are "hard-coded" (except P) by using the HSP lab
*  data.
*
*  FORTRAN Name: VPHAIN.FOR
*
*  Keywords of Accessed Files and Tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       None
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       CDINTR
*  SDAS:
*       UMSPUT
*  OTHERS:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       04-15-86     J.-C. Hsu        Design and coding
*     2       07-20-86     J.-C. Hsu        change the algorithm of calculating
*                                           P, Q, and R
*     3       08-20-87     J.-C. Hsu        F77 standard
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --number of input data points
        INTEGER                 NPTS
*                                       --threshold setting 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:
*                                       --index of highest cout rate
        INTEGER                 IYMAX,
     :                          STATOK,
*                                       --loop indices
     :                          I
*                                       --highest count rate
        DOUBLE PRECISION        YMAX,
     :                          X1, X2
*                                       --error message context
        CHARACTER*130           CONTXT
*=========================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 = 6
        NVAR = 2
*
        COEFF(4) = 1000.D0
        COEFF(5) = 3.D0
        COEFF(6) = 0.D0
*
*  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
*
*  find the highest count rate, and make sure it's at a small discriminator
*  threshold setting (<5)
*
        YMAX = 0
        IYMAX = 0
        DO 10 I = 1, NPTS
            IF (Y(I) .GT. YMAX) THEN
                YMAX = Y(I)
                IYMAX = I
            END IF
   10   CONTINUE
*
        IF (X(IYMAX) .GT. 4) THEN
            STATUS = ERRNUM(1)
            CONTXT = 'lack of data points near the origin'
            GO TO 999
        END IF
*
*  find Q (=COEFF(2)), which is approximately at half total area
*
        CALL CDINTR (X, Y, NPTS, YMAX/2.,
     :                  COEFF(2))
        COEFF(1) = YMAX / 2.506628D0
*
        CALL CDINTR (X, Y, NPTS, YMAX * 0.15865D0,
     :                  X1)
        CALL CDINTR (X, Y, NPTS, YMAX * 0.84135D0,
     :                  X2)
        COEFF(3) = ABS(X2 - X1) / 2.D0
*
        STATUS = OK
        GO TO 1000
*
*  write error message
*
  999   CALL UMSPUT ('VPHAIN: ' // CONTXT, DEST, PRIO,   STATOK)
*
 1000   RETURN
        END