SUBROUTINE CDCTRD (
*
*  inputs
*
     :                  COUNT, FLAG, X, Y, OKVAL, HITVAL, DIM, NPEAK,
     :                  HPT, VPT, HSTPT, VSTPT,
*
*  output
*
     :                  PEAKCT, XCENT, YCENT, DXCENT, DYCENT, STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  Find the brightest NPEAK points and calculate their centroid coordinates
*
*  Description:
*  ------------
*  Maximum size of the array is 50 by 50.
*
*  FORTRAN name: CDCTRD.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       01-20-87     J.-C. HSU        design and coding
*     2       11-05-87     J.-C. HSU        F77 SDAS
*-------------------------------------------------------------------------------
*
*== input:
*                                       --number of columns and rows of the
*                                       --area scan
        INTEGER         HPT, VPT,
*                                       --number of deflection steps in horizontal
*                                       --and vertical directions
     :                  HSTPT, VSTPT,
*                                       --dimension of COUNT, FLAG, X, and Y
*                                       --as declared in the calling routine
     :                  DIM,
*                                       --number of highest counts to be picked
*                                       --out
     :                  NPEAK
*                                       --observed (average) number of counts,
*                                       --and its flag
        REAL            COUNT(DIM, DIM), FLAG(DIM, DIM),
*                                       --coordinates of the area scan points
     :                  X(DIM, DIM), Y(DIM, DIM),
*                                       --mask value of good data points
     :                  OKVAL,
*                                       --mask value of particle events
     :                  HITVAL
*
*== output
*                                       --highest count
        REAL            PEAKCT,
*                                       --centroid coordinates of the highest
*                                       --NPEAK points and its standard
*                                       --deviation
     :                  XCENT, YCENT, DXCENT, DYCENT
*                                       --error status
        INTEGER         STATUS
*
*== local:
*                                       --maximum number and sum
        REAL            MAX, SUM, XSUM, YSUM, XSUMSQ, YSUMSQ,
*                                       --local flag array
     :                  XFLAG(50, 50)
*                                       --error number and loop indices
        INTEGER         STATOK, I, J, K, JJ, KK
        CHARACTER*5     CHAR5
*                                       --error message context
        CHARACTER*130   CONTXT, MESS
*=========================begin hsp.inc=========================================
*                                       --status return code
        INTEGER         OK, ERRNUM(20)
        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/
*                                       --message destination and priority
        DATA DEST, PRIO /1, 0/
*=========================end hsp.inc===========================================
*------------------------------------------------------------------------------
*
*  check the dimension of XFLAG
*
        IF (DIM .GT. 50) THEN
            STATUS = ERRNUM(1)
            CONTXT = 'improper local array dimension declaration'
            GO TO 999
        END IF
*
*  equate the local flag to the input flag
*
        DO 20 I = 1, DIM
            DO 10 J = 1, DIM
                XFLAG(I, J) = FLAG(I, J)
   10       CONTINUE
   20   CONTINUE
*
*  look for the NPEAK highest points
*
        SUM = 0.
        XSUM = 0.
        YSUM = 0.
        XSUMSQ = 0.
        YSUMSQ = 0.
*
        DO 50 I = 1, NPEAK
            MAX = -1.
*
*  process each row
*
            DO 40 J = 1, VPT
*
*  process each column
*
                DO 30 K = 1, HPT
                    IF (XFLAG(J, K) .NE. HITVAL .AND.
     :                          COUNT(J, K) .GT. MAX) THEN
                        MAX = COUNT(J, K)
                        JJ = J
                        KK = K
                    END IF
   30           CONTINUE
   40       CONTINUE
*
*  if no maximum number
*
            IF (MAX .LT. 0) THEN
                STATUS = ERRNUM(2)
                WRITE (CHAR5, '(I5)') NPEAK
                CONTXT = 'good data points is less than ' // CHAR5
                GO TO 999
            END IF
            XFLAG(JJ, KK) = HITVAL
*
            SUM = MAX + SUM
            XSUM = MAX * X(JJ, KK) + XSUM
            YSUM = MAX * Y(JJ, KK) + YSUM
            XSUMSQ = MAX * X(JJ, KK) ** 2 + XSUMSQ
            YSUMSQ = MAX * Y(JJ, KK) ** 2 + YSUMSQ
*
*  pick out the highest count
*
            IF (I .EQ. 1) THEN
                PEAKCT = MAX
            END IF
   50   CONTINUE
*
*  centroid coordinate and its standard deviation
*
        XCENT = XSUM / SUM
        YCENT = YSUM / SUM
        DXCENT = FLOAT(HSTPT)
        DYCENT = FLOAT(VSTPT)
*
        STATUS = OK
        GO TO 1000
*
*  write error message
*
  999   MESS = 'CDCTRD: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,   STATOK)
*
 1000   RETURN
        END