SUBROUTINE CDROTM (
*
*  inputs
*
     :                  MJD,
*
*  outputs
*
     :                  MATRIX)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  compute the precession matrix from J2000.0 TO MJD
*
*  Description:
*  ------------
*  See Astronomical Almanac 1984, p. S19 for equations used here.
*  Adapted from IRAF noao ast_rotmatrix
*
*  FORTRAN name: CDROTM.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_rotmatrix
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --epochs in Modified Julian days
        DOUBLE PRECISION        MJD
*
*== output:
*                                       --precession rotation matrix
        DOUBLE PRECISION        MATRIX(3, 3)
*
*== local:
*                                       --Julian day of epoch J2000.0
        DOUBLE PRECISION        MJ2000,
*                                       --centuries since J2000
     :                          T,
*                                       --degrees per radian
     :                          RADN,
*                                       --rotation angles and their sines and
*                                       --cosines
     :                          A, B, C, SA, SB, SC, CA, CB, CC
        PARAMETER               (MJ2000 = 51544.5D0)
*
*------------------------------------------------------------------------------
*
*  degrees of one radian
*
        RADN = 45.D0 / ATAN(1.D0)
*
*  number of centuries from J2000
*
        T = (MJD - MJ2000) / 36525.D0
*
*  calculate zeta(A), zee(A), theta(A) in seconds of arc
*
        A = T * (2306.2181D0 + T * ( 0.30188D0 + T * 0.017998))
        B = T * (2306.2181D0 + T * ( 1.09468D0 + T * 0.018203))
        C = T * (2004.3109D0 + T * (-0.42665D0 - T * 0.041833))
*
*  convert angles to radians
*
        A = A / (3600.D0 * RADN)
        B = B / (3600.D0 * RADN)
        C = C / (3600.D0 * RADN)
*
*  calculate sines and cosines of angles
*
        CA = COS (A)
        SA = SIN (A)
        CB = COS (B)
        SB = SIN (B)
        CC = COS (C)
        SC = SIN (C)
*
*  compute the rotation matrix from sines and cosines of the rotation angles
*
        MATRIX(1,1) = CA * CB * CC - SA * SB
        MATRIX(2,1) = -SA * CB * CC -CA * SB
        MATRIX(3,1) = -CB * SC
        MATRIX(1,2) = CA * SB * CC + SA * CB
        MATRIX(2,2) = -SA * SB * CC + CA * CB
        MATRIX(3,2) = - SB * SC
        MATRIX(1,3) = CA * SC
        MATRIX(2,3) = - SA * SC
        MATRIX(3,3) = CC
*
        RETURN
        END