        SUBROUTINE CDPRC (
*
*  inputs
*
     :                  XYZ1, MJD1, MJD2,
*
*  outputs
*
     :                  XYZ2)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  precess rectangular coordinates from MJD1 to MJD2
*
*  Description:
*  ------------
*  See Astronomical Almanac 1984, p. S19 for equations used here.
*  Adapted from IRAF noao ast_precess
*
*  FORTRAN name: CDPRC.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       None
*  SDAS:
*       None
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       02-28-88     J.-C. HSU        coding, adapted from noao ast_mst
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --rectangular coordinates at epoch MJD1
        DOUBLE PRECISION        XYZ1(3),
*                                       --epochs in Modified Julian days
     :                          MJD1, MJD2
*
*== output:
*                                       --rectangular coordinates at epoch MJD2
        DOUBLE PRECISION        XYZ2(3)
*
*== local:
*                                       --Modified Julian day of epoch J2000.0
        DOUBLE PRECISION        MJ2000,
*
     :                          XYZ(3),
*                                       --rotation matrix
     :                          P(3, 3)
        INTEGER                 I
        PARAMETER               (MJ2000 = 51544.5D0)
*
*------------------------------------------------------------------------------
*
        IF (MJD1 .EQ. MJD2) THEN
            DO 10 I = 1, 3
                XYZ2(I) = XYZ1(I)
   10       CONTINUE
        ELSE
*
*  if MJD1 is not the standard epoch (J2000.), precess from MJD1 to J2000
*
            IF (MJD1 .EQ. MJ2000) THEN
                DO 20 I = 1,3
                    XYZ(I) = XYZ1(I)
   20           CONTINUE
            ELSE
*
*  precess from MJD1 TO J2000, so the rotation matrix is the inverse of P
*  which is the transpose of P
*
                CALL CDROTM (MJD1,    P)
*
                DO 30 I = 1, 3
                    XYZ(I) = P(I,1) * XYZ1(1) + P(I,2) * XYZ1(2) +
     :                       P(I,3) * XYZ1(3)
   30           CONTINUE
            END IF
*
*  if MJD2 is not the standard epoch (J2000.), precess from J2000 to MJD2
*
            IF (MJD2 .EQ. MJ2000) THEN
                DO 40 I = 1,3
                    XYZ2(I) = XYZ(I)
   40           CONTINUE
            ELSE
*
                CALL CDROTM (MJD2,    P)
*
                DO 50 I = 1, 3
                    XYZ2(I) = P(1,I) * XYZ(1) + P(2,I) * XYZ(2) +
     :                        P(3,I) * XYZ(3)
   50           CONTINUE
            END IF
        END IF
*
        RETURN
        END
