SUBROUTINE VSTATE (
*
*  inputs
*
     :                          T, T0, A, P,
*
*  outputs
*
     :                          X, DXDT)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  calculate the state vector and their time derivatives
*
*  Description:
*  ------------
*  Using the formula described in page B-25 of SOGS SE-06-2, 19 December 1984.
*
*  FORTRAN name: VSTATE.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       05-28-89     J.-C. HSU        coding
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --epoch and time zero point in seconds
*                                       --from 1980.0
        DOUBLE PRECISION        T, T0,
*                                       --orbit coefficients
     :                          A(3), P(4, 7)
*
*== output:
*                                       --state vector of the desired epoch
*                                       --and its time derivatives (in km
*                                       --and km/sec)
        DOUBLE PRECISION        X(3), DXDT(3)
*
*== local:
*
        INTEGER                 J
        DOUBLE PRECISION        X1, X2, X3, DX1, DX2, DX3, DT, G,
     :                          DGDT, F(4), DF(4), SING, COSG, SIN2G,
     :                          COS2G, SINF1, COSF1, SINF3, COSF3,
     :                          SINF4, COSF4, RADN
        PARAMETER               (RADN = 57.29577951308232 D0)
*
*------------------------------------------------------------------------------
*=========================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===========================================
*
*  calculate g, and its time derivative
*
        DT = T - T0
*
        G = A(1) + (A(2) + (A(3) * DT)) * DT
        DGDT = A(2) + 2.D0 * A(3) * DT
*
*  convert degrees to radians
*
        G = G / RADN
        DGDT = DGDT / RADN
*
        SING = SIN(G)
        COSG = COS(G)
        SIN2G = 2.D0 * SING * COSG
        COS2G = COSG * COSG - SING * SING
*
*  calculate Fj's and their time derivatives
*
        DO 10 J = 1, 4
            F(J) = P(J,1) + (P(J,2) + (P(J,3) * DT)) * DT +
     :             P(J,4) * COSG + P(J,5) * SING +
     :             P(J,6) * COS2G + P(J,7) * SIN2G
*
            DF(J) = P(J,2) + 2.D0 * P(J,3) * DT +
     :                DGDT * (P(J,5) * COSG - P(J,4) * SING +
     :                2.D0 * P(J,7) * COS2G - 2.D0 * P(J,6) * SIN2G)
   10   CONTINUE
*
*  convert F(1), F(3), F(4) from degrees to radians and add g for F1
*
        F(1) = F(1)/RADN + G
        F(3) = F(3)/RADN 
        F(4) = F(4)/RADN 
        DF(1) = DF(1)/RADN + DGDT
        DF(3) = DF(3)/RADN 
        DF(4) = DF(4)/RADN 
*
*  evaluate the state vector and its time derivative
*
        SINF1 = SIN(F(1))
        COSF1 = COS(F(1))
        SINF3 = SIN(F(3))
        COSF3 = COS(F(3))
        SINF4 = SIN(F(4))
        COSF4 = COS(F(4))
*
        X1 = COSF3 * COSF1 - SINF3 * COSF4 * SINF1
        X2 = SINF3 * COSF1 + COSF3 * COSF4 * SINF1
        X3 = SINF4 * SINF1
*
        DX1 = -SINF3 * DF(3) * COSF1 - COSF3 * SINF1 * DF(1) -
     :         COSF3 * DF(3) * COSF4 * SINF1 +
     :         SINF3 * SINF4 * DF(4) * SINF1 -
     :         SINF3 * COSF4 * COSF1 * DF(1)
        DX2 =  COSF3 * DF(3) * COSF1 - SINF3 * SINF1 * DF(1) -
     :         SINF3 * DF(3) * COSF4 * SINF1 -
     :         COSF3 * SINF4 * DF(4) * SINF1 +
     :         COSF3 * COSF4 * COSF1 * DF(1)
        DX3 =  COSF4 * DF(4) * SINF1 + SINF4 * COSF1 * DF(1)
*
        X(1) = F(2) * X1
        X(2) = F(2) * X2
        X(3) = F(2) * X3
*
        DXDT(1) = DF(2) * X1 + F(2) * DX1
        DXDT(2) = DF(2) * X2 + F(2) * DX2
        DXDT(3) = DF(2) * X3 + F(2) * DX3
*
        RETURN
        END