      SUBROUTINE TRUK2VEL(T,A,V1,V2,DVELDA)
C----------------------------------------------------------------
C THIS SUBROUTINE CALCULATES RADIAL VELOCITIES  
C FOR A GIVEN TIME AND A GIVEN SET OF ORBITAL ELEMENTS.  
C THE FORMULAS ARE BASED ON GREEN, SPHERICAL ASTRONOMY, CH. 19.
C NOTE THAT, UNLIKE MOST OF THE SUBROUTINES IN THIS PACKAGE, THIS
C ONE DOES  N O T  USE THE INCLUDE FILES 'PARAMETR.INC' AND 
C 'ELLIPSE.INC'.
C
C     23 August 1991      J.T. Armstrong
C      1 July   1994      J.T. Armstrong:  modified from TRUE2VEL to
C                             use K1, K2 instead of a1 sin i, a2 sin i
C
C----------------------------------------------------------------
C THE BASIC FORMULA FOR STAR i IS:
C
C     V(rad) = V0 + Ki * [ cos( V(t) + omega ) + e cos(omega) ]
C                 + dV0/dt * (t - T)
C
C with
C
C     Ki = (semi-)amplitude of velocity curve
C
C     V(t) = true anomaly
C
C     omega = argument of periastron
C
C     e = eccentricity
C
C     T = epoch of periastron
C
C     V0 = mean radial velocity
C
C     dV0/dt = acceleration (km/sec/yr)
C
C Note that for star 1, use omega, but for star 2, use omega+pi.

C VARIABLES:
C     M_ANOM      MEAN ANOMALY
C     E_ANOM      ECCENTRIC ANOMALY
C     T_ANOM      TRUE ANOMALY
C
C     A           ARRAY CONTAINING ORBITAL ELEMENTS
C     A_TRUE      SEMIMAJOR AXIS OF TRUE RELATIVE ORBIT (mas); NOT USED
C     K1, K2      VELOCITY AMPLITUDES FOR STARS 1, 2
C     E_TRUE      ECCENTRICITY OF TRUE ORBIT
C     PERIOD      DAYS
C     EPOCH       JULIAN DATE
C     INCLIN      INCLINATION ANGLE (rad); NOT USED, SUBSUMED IN A1,2_SINI
C     ASC_NODE    ANGLE OF ASCENDING NODE (rad); NOT USED
C     ARG_PERI    ARGUMENT OF PERIASTRON (rad)
C     V0          RADIAL VELOCITY OF CENTER OF MASS (km/s)
C     DV0DT       ACCELERATION OF RADIAL VELOCITY OF C.M. (km/s/yr)
C
C     T           OBSERVATION TIME (Julian days)
C
C     A1_SINI     SEMIMAJOR AXIS OF TRUE ORBIT OF STAR 1 * SIN(I) (km/s)
C     A2_SINI     SEMIMAJOR AXIS OF TRUE ORBIT OF STAR 2 * SIN(I) (km/s)
C     V1, V2      CALCULATED RADIAL VELOCITIES OF STAR 1, 2 (km/s)
C
C     DV1DA, DV2DA:     DERIVATIVES OF V1, V2 W.R.T. ELEMENTS A
C     DV1DK1, DV1DK2:   DERIVATIVES OF V1, V2 W.R.T. K1, K2
C     DEDE, DEDP, DEDT: DERIVATIVES OF ECCENTRIC ANOMALY
C     DVDE, DVDP, DVDT: DERIVATIVES OF TRUE ANOMALY
C
C     VELC        SPEED OF LIGHT (m/s)
C     SECPDAY     SECONDS PER DAY
C
C     J           DO LOOP INDEX
C     SUM         USED AS DUMMY
C

      IMPLICIT UNDEFINED (A-Z)

      REAL*8    M_ANOM, E_ANOM, T_ANOM, SUM
      REAL*8    A_TRUE,E_TRUE,PERIOD,EPOCH,INCLIN,ASC_NODE,ARG_PERI
      REAL*8    A(11), T, PI, TWOPI, RPDEG, VELC, SECPDAY, JYEAR
      REAL*8    A1_SINI,A2_SINI,V0,DV0DT

      REAL*8    K1,K2,V1,V2,DVELDA(11,2)
      REAL*8    COSV2, SINV_AP
      REAL*8    E_COSE, SEC_SEC
      REAL*8    DEDE, DEDP, DEDT, DVDE, DVDP, DVDT
      REAL*8    DV1DV0,DV1DK1,DV1DP,DV1DE,DV1DAP,DV1DT,DV1DACC
      REAL*8    DV2DV0,DV2DK2,DV2DP,DV2DE,DV2DAP,DV2DT,DV2DACC

      REAL*8    TMOD

      INTEGER*4 J

      PARAMETER (PI = 3.14159265359D0)
      PARAMETER (TWOPI = 2.0D0 * PI)
      PARAMETER (RPDEG = PI/180D0)
      PARAMETER (VELC = 2.99792458D8)
      PARAMETER (SECPDAY = 60D0 * 60D0 *24D0)
      PARAMETER (JYEAR = 365.25)

C      TWOPI = 2.0D0 * PI

C     WRITE(6,'(A)') ' TRUE2VEL:  ELEMENTS TO TRUE2VEL:'
C     WRITE(6,
C           '(A,F9.5,A,F8.5,A,F8.3,A,F12.4 / A,F10.5,A,F9.4,A,F9.4)')
C    1      '  A: ',A(1),'  E: ',A(2),'  P: ',A(3),'  T: ',A(4),
C    2      '  I: ',A(5)/RPDEG,' ASC: ',A(6)/RPDEG,'  ARG: ',
C    3      A(7)/RPDEG
C     WRITE(6,'(3(A,F9.5))') ' K1: ',A(8),'  K2: ',
C    1      A(9),'  V0: ',A(10),'  dV0/dt:  ',A(11)

      A_TRUE      = A(1)
      E_TRUE      = A(2)
      PERIOD      = A(3)
      EPOCH       = A(4)
      INCLIN      = A(5)
      ASC_NODE    = A(6)
      ARG_PERI    = A(7)
      K1          = A(8)
      K2          = A(9)
      V0          = A(10)
      DV0DT       = A(11)

C THE MEAN ANOMALY IN RADIANS IS

D     TMOD        = DMOD(T-EPOCH,PERIOD)
D     WRITE(6,*) ' TMOD = ',TMOD
      M_ANOM      = TWOPI* (DMOD(T-EPOCH, PERIOD) / PERIOD)

C GET THE FIRST ESTIMATE OF THE ECCENTRIC ANOMALY

      E_ANOM      = M_ANOM + E_TRUE * SIN(M_ANOM)
     1                   + 0.5*E_TRUE**2 * SIN(2.*M_ANOM)

C ITERATE TO GET THE ECCENTRIC ANOMALY

      DO J = 1, 5
          E_ANOM  = E_ANOM + (M_ANOM-E_ANOM+E_TRUE*SIN(E_ANOM))/
     1                         (1. - E_TRUE * COS(E_ANOM))
      END DO

C GET THE TRUE ANOMALY

      T_ANOM      = 2*ATAN(
     1                 SQRT((1.+E_TRUE)/(1.-E_TRUE))
     2                 * TAN(E_ANOM/2.)
     3                )

C     WRITE(6,'(1X,A,F7.2,2X,F7.2,2X,F7.2)') 
C    1         'TRUE2VEL: MEAN, ECC, TRUE ANOM.S BEFORE CHECK: ',
C    2         M_ANOM,E_ANOM,T_ANOM

C SINCE E_ANOM IS ALWAYS POSITIVE, T_ANOM SHOULD BE TOO:

      IF (T_ANOM.LT.0.) THEN
          T_ANOM = TWOPI+DMOD(T_ANOM,TWOPI)
      ELSE
          T_ANOM = DMOD(T_ANOM,TWOPI)
      END IF

C     WRITE(6,'(1X,A,F11.3,2X,F7.2,2X,F7.2,2X,F7.2)') 
C    1         'TRUE2VEL: T AND MEAN, ECC, TRUE ANOM.S:  ',
C    2         T,M_ANOM,E_ANOM,T_ANOM

C CALCULATE A1 SIN I, A2 SIN I:

      A1_SINI = K1 * SQRT(1-E_TRUE**2) * PERIOD * SECPDAY / TWOPI

      A2_SINI = K2 * SQRT(1-E_TRUE**2) * PERIOD * SECPDAY / TWOPI

C CALCULATE THE RADIAL VELOCITIES:
C Note that using the minus sign in the expression for V2 corresponds to
C adding PI to ARG_PERI.

      V1           = V0 + K1 * 
     1                  ( COS(T_ANOM+ARG_PERI) + E_TRUE*COS(ARG_PERI) )
     2                  + DV0DT * (T - EPOCH)/JYEAR   
      V2           = V0 - K2 * 
     1                  ( COS(T_ANOM+ARG_PERI) + E_TRUE*COS(ARG_PERI) )
     2                  + DV0DT * (T - EPOCH)/JYEAR   


C EVALUATE THE DERIVATIVES
C FIRST, SOME USEFUL EXPRESSIONS

      COSV2        = COS(T_ANOM/2.)
      SINV_AP      = SIN(T_ANOM+ARG_PERI)

      E_COSE       = 1. - E_TRUE * COS(E_ANOM)
      SEC_SEC      = (COSV2/COS(E_ANOM/2.))**2
     1                    * SQRT( (1.+E_TRUE) / (1.-E_TRUE) )

C DERIVATIVES OF THE ECCENTRIC ANOMALY E_ANOM W.R.T. ECCENTRICITY,
C PERIOD, AND EPOCH:

      DEDE         = SIN(E_ANOM) / E_COSE

      DEDP         = -TWOPI * (T - EPOCH) / 
     1                      (PERIOD**2 * E_COSE)

      DEDT         = -TWOPI / (PERIOD * E_COSE)

C DERIVATIVE OF THE TRUE ANOMALY T_ANOM W.R.T. ECCENTRICITY, PERIOD,
C AND EPOCH:

      DVDE         = SEC_SEC * DEDE +
     1                   2. * TAN(E_ANOM/2.) * COSV2**2 *
     2                   SQRT( (1.-E_TRUE)/(1.+E_TRUE) ) / 
     3                   (1.-E_TRUE)**2

      DVDP         = SEC_SEC * DEDP

      DVDT         = SEC_SEC * DEDT

C DERIVATIVES OF V1 W.R.T. V0, K1, PERIOD, ECCENTRICITY, ARG. OF
C PERIASTRON, EPOCH, AND ACCELERATION:

      DV1DV0       = 1.0D0

      DV1DK1       = (V1 - V0) / K1

      DV1DP        = -(V1 - V0)/PERIOD - K1*SINV_AP*DVDP

      DV1DE        = ( E_TRUE*(V1-V0) ) / SQRT(1.-E_TRUE**2) +
     1                 K1 * ( COS(ARG_PERI) - SINV_AP*DVDE )

      DV1DAP       = -K1 * ( SINV_AP + E_TRUE*SIN(ARG_PERI) )

      DV1DT        = -K1 * SINV_AP * DVDT

      DV1DACC      = (T - EPOCH)/JYEAR

C DERIVATIVES OF V2 W.R.T. V0, K2, PERIOD, ECCENTRICITY, ARG. OF
C PERIASTRON, EPOCH, AND ACCELERATION:
C Note that using minus/plus signs here in place of plus/minus signs in
C the expressions for the derivatives of V1 correspond to adding PI to 
C ARG_PERI.

      DV2DV0       = 1.0D0

      DV2DK2       = (V2 - V0) / K2

      DV2DP        = -(V2 - V0)/PERIOD + K2*SINV_AP*DVDP

      DV2DE        = ( E_TRUE*(V2-V0) ) / SQRT(1.-E_TRUE**2) -
     1                 K2 * ( COS(ARG_PERI) - SINV_AP*DVDE )

      DV2DAP       = K2 * ( SINV_AP + E_TRUE*SIN(ARG_PERI) )

      DV2DT        = K2 * SINV_AP * DVDT

      DV2DACC      = (T - EPOCH)/JYEAR

C NOW PUT THE DERIVATIVES INTO THE DVELDA ARRAY:

      DVELDA(1,1)      = 0.
      DVELDA(2,1)      = DV1DE
      DVELDA(3,1)      = DV1DP
      DVELDA(4,1)      = DV1DT
      DVELDA(5,1)      = 0.
      DVELDA(6,1)      = 0.
      DVELDA(7,1)      = DV1DAP
      DVELDA(8,1)      = DV1DK1
      DVELDA(9,1)      = 0.
      DVELDA(10,1)     = DV1DV0
      DVELDA(11,1)     = DV1DACC
      
      DVELDA(1,2)      = 0.
      DVELDA(2,2)      = DV2DE
      DVELDA(3,2)      = DV2DP
      DVELDA(4,2)      = DV2DT
      DVELDA(5,2)      = 0.
      DVELDA(6,2)      = 0.
      DVELDA(7,2)      = DV2DAP
      DVELDA(8,2)      = 0.
      DVELDA(9,2)      = DV2DK2
      DVELDA(10,2)     = DV2DV0
      DVELDA(11,2)     = DV2DACC

C     WRITE(6,'(A)') ' TRUE2VEL:'
C     WRITE(6,'(A)')
C    1     '       epoch       mean anom.   ecc. anom.   true anom.'
C RULER LINE:
CJD ttttttttttttXXXmmmmmmmmmmXXXeeeeeeeeeeXXXvvvvvvvvvv
C     WRITE(6,'(1X,A,F12.3,3(3X,F10.5))')
C    1    'JD ',T-2400000., M_ANOM*PERIOD/TWOPI, E_ANOM*PERIOD/TWOPI,
C    2    T_ANOM*PERIOD/TWOPI

C     WRITE(6,'(A,F7.3,2X,F7.3)') ' A1_SINI, A2_SINI = ',
C    1    A1_SINI, A2_SINI
C     WRITE(6,'(A)')

C     WRITE(6,'(1X,A,A)')
C    1    '        dv/de       dv/dP       dv/dT',
C    2    '     dv/domega  dv/d(a sini)   dv/dv0'
C RULER LINE:
CX...........X...........X...........X...........X...........X...........
C     WRITE(6,'(1X,A,6(1X,G11.3))') 'V1: ',
C    1    DV1DE,DV1DP,DV1DT,DV1DAP,DV1DK1,DV1DV0
C     WRITE(6,'(1X,A,6(1X,G11.3))') 'V2: ',
C    1    DV2DE,DV2DP,DV2DT,DV2DAP,DV2DK2,DV2DV0

      
      RETURN

      END

