        SUBROUTINE CDCIRC (
*
*  inputs:
*
     :                          X, Y, DX, DY, NPTS, SCHEME, FRAC, DIM,
*
*  inputs/outputs:
*
     :                          ITER, TOLERN,
*
*  outputs:
*
     :                          NCOEFF, NVAR, COEFF, MATRIX, CHISQ,
     :                          STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  Fit a circle to a set of contour points.
*
*  Description:
*  ------------
*  Maximum number of input data points is 2000.
*
*  FORTRAN name: CDCIRC.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       CDCLIN, CDCLLS
*  SDAS:
*       UMSPUT
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       01-20-87     J.-C. HSU        design and coding
*     2       05-20-87     J.-C. HSU        F77 SDAS
*-------------------------------------------------------------------------------
*
*== inputs:
*                                       --X and Y arrays and their standard
*                                       --deviations
        REAL                    X(1), Y(1), DX(1), DY(1),
*                                       --weighting schemes
     :                          SCHEME(1),
*                                       --specified fraction applied to the
*                                       --coefficient modification
     :                          FRAC
*                                       --total number of valid input data points
        INTEGER                 NPTS,
*                                       --dimension of MATRIX
     :                          DIM
*
*== inputs/outputs:
*
*                                       --maximum number of iterations of least
*                                       --square fitting
        INTEGER                 ITER
*                                       --tolerance of sigma-squared difference
        REAL                    TOLERN
*
*== outputs:
*                                       --number of coefficients and variables
        INTEGER                 NCOEFF, NVAR,
*                                       --error status
     :                          STATUS
*                                       --fitting coefficients, covariance
*                                       --matrix, and chi-squared
        DOUBLE PRECISION        COEFF(*), MATRIX(DIM, DIM), CHISQ
*
*== local:
*                                       --sum of X and Y
        REAL                    XSUM, YSUM
*                                       --X and Y arrays and their standard
*                                       --deviations in double precision
        DOUBLE PRECISION        XX(2000), YY(2000), DXX(2000),
     :                          DYY(2000),
*                                       --X, Y offset
     :                          XSHIFT, YSHIFT
*                                       --loop index
        INTEGER                 I,
*                                       --error status
     :                          STATOK
*                                       --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===========================================
*------------------------------------------------------------------------------
*
        IF (NPTS .GT. 2000) THEN
            STATUS = ERRNUM(1)
            CONTXT = 'more than 2000 input data points'
            GO TO 999
        END IF
*
*  get initial values of the coefficients
*
        CALL CDCLIN (X, Y, NPTS,    COEFF, NCOEFF, NVAR, STATUS)
        IF (STATUS .NE. OK) THEN
            CONTXT = 'cannot get initial coefficients'
            GO TO 999
        END IF
*
*  shift the origin such taht the mean coordinates are zero
*
        XSUM = 0.
        YSUM = 0.
*
        DO 10 I = 1, NPTS
            XSUM = X(I) + XSUM
            YSUM = Y(I) + YSUM
   10   CONTINUE
*
        XSHIFT = DBLE(XSUM / REAL(NPTS))
        YSHIFT = DBLE(YSUM / REAL(NPTS))
*
*  convert input data to double precision
*
        DO 20 I = 1, NPTS
            XX(I) = DBLE(X(I)) - XSHIFT
            YY(I) = DBLE(Y(I)) - YSHIFT
            DXX(I) = DBLE(DX(I))
            DYY(I) = DBLE(DY(I))
   20   CONTINUE
*
*  do the least square fitting
*
        CALL CDCLLS (XX, YY, DXX, DYY, NPTS, SCHEME, FRAC,
     :                  NCOEFF, DIM, NVAR,
     :                  COEFF, TOLERN, ITER,
     :                  CHISQ, MATRIX, STATUS)
        IF (STATUS .NE. OK) THEN
            CONTXT = 'error in circle least-square fitting'
            GO TO 999
        END IF
*
*  shift the center back to the original coordinate system
*
        COEFF(1) = COEFF(1) + XSHIFT
        COEFF(2) = COEFF(2) + YSHIFT
*
        STATUS = OK
        GO TO 1000
*
*  write error message
*
  999   MESS = 'CDCIRC: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,   STATOK)
*
 1000   RETURN
        END
