      SUBROUTINE TRUE2APP(T,A,REST,THEST,RAEST,DECEST,DYDA,PR)
C----------------------------------------------------------------
C THIS SUBROUTINE CALCULATES THE SEPARATIONS AND POSITION ANGLES 
C FOR GIVEN TIMES AND FOR A GIVEN SET OF ORBITAL ELEMENTS.  
C THE FORMULAS ARE BASED ON W.D. HEINTZ, DOUBLE STARS, AROUND P. 24.
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 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 ORBIT
C     E_TRUE      ECCENTRICITY OF TRUE ORBIT
C     PERIOD
C     EPOCH
C     INCLIN      INCLINATION ANGLE
C     ASC_NODE    ANGLE OF ASCENDING NODE
C     ARG_PERI    ARGUMENT OF PERIASTRON
C
C     T           OBSERVATION TIME IN JULIAN DAYS
C     R_TRUE      DIST. BTW. STARS IN TRUE ORBIT
C     THEST       POSITION ANGLE IN APPARENT ORBIT
C     REST        SEPARATION IN APPARENT ORBIT
C     RAEST       R.A. SEPARATION IN APPARENT ORBIT
C     DECEST      DEC. SEPARATION IN APPARENT ORBIT
C
C     DYDA        DERIVATIVES OF REST, THEST W.R.T. ELEMENTS A
C     DEDE, DEDP, DEDT: DERIVATIVES OF ECCENTRIC ANOMALY
C     DVDE, DVDP, DVDT: DERIVATIVES OF TRUE ANOMALY
C     DRDA, DRDE, DRDP, DRDT: DERIVATIVES OF TRUE SEPARATION
C     DTHETADA, ETC:    DERIVATIVES OF APPARENT POSITION ANGLE
C     DRHODA, ETC:      DERIVATIVES OF APPARENT SEPARATION
C
C     J           DO LOOP INDEX
C     SUM         USED AS DUMMY
C
C     PR          UNIT NUMBER OF OUTPUT FILE 

      IMPLICIT UNDEFINED (A-Z)

      REAL*8    M_ANOM, E_ANOM, T_ANOM, R_TRUE, SUM
      REAL*8    M_TEST1, M_TEST2, E_TEST
      REAL*8    A_TRUE,E_TRUE,PERIOD,EPOCH,INCLIN,ASC_NODE,ARG_PERI
      REAL*8    A(7), DYDA(7,2), T, PI, TWOPI, RPDEG

      REAL*8    THEST,REST,RAEST,DECEST
      REAL*8    DV_DTIME,DTH_DTIME,DR_DT1,DR_DT2,DR_DTIME
      REAL*8    DTH_DHR,DR_DHR
      REAL*8    COSV2, COSV_AP, SINV_AP, COSTH_AN, COSI, SINI, SINV 
      REAL*8    TANV_AP, E_COSE, E_COSV, SEC_SEC, R_RHO, R2_RHO
      REAL*8    COS2_COS2SIN2, COSSINSIN2, SIN_COS
      REAL*8    DEDE, DEDP, DEDT, DVDE, DVDP, DVDT
      REAL*8    DRDA, DRDE, DRDP, DRDT
      REAL*8    DTHETADA,DTHETADE,DTHETADP,DTHETADT,DTHETADI
      REAL*8    DTHETADAN, DTHETADAP
      REAL*8    DRHODA,DRHODE,DRHODI,DRHODP,DRHODT,DRHODAN,DRHODAP

      INTEGER*4 J,PR

      PARAMETER (PI = 3.14159265359D0)
      PARAMETER (RPDEG = PI/180D0)

      TWOPI = 2.0D0 * PI

C     WRITE(PR,'(A)') ' TRUE2APP:  ELEMENTS TO TRUE2APP:'
C     WRITE(PR,'(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

      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)

C THE MEAN ANOMALY IN RADIANS IS

      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 TEST TO SEE IF WE'VE GOTTEN E_ANOM RIGHT:

D     M_TEST1 = E_ANOM - E_TRUE*SIN(E_ANOM)

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(PR,'(1X,A,F7.2,2X,F7.2,2X,F7.2)') 
C    1         'TRUE2APP: 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 TEST TO SEE IF WE'VE GOTTEN T_ANOM RIGHT:

D     E_TEST = 2.*ATAN(SQRT((1.-E_TRUE)/(1.+E_TRUE))*TAN(T_ANOM/2.))
D     M_TEST2 = E_TEST - E_TRUE*SIN(E_TEST)

C     WRITE(PR,'(1X,A,F11.3,5(2X,F7.2))') 
C    1         'TRUE2APP: T AND MEAN, ECC, TRUE ANOM.S:  ',
C    2         T,M_ANOM/RPDEG,E_ANOM/RPDEG,T_ANOM/RPDEG,
C    3         M_TEST1/RPDEG,M_TEST2/RPDEG,E_TEST/RPDEG

C GET THE DISTANCE BETWEEN STARS (STILL IN THE TRUE ORBIT)

      R_TRUE      = A_TRUE * (1. - E_TRUE**2) / 
     1                     (1. + E_TRUE*COS(T_ANOM))

C      WRITE (PR,*) 'TRUE2APP: R_TRUE = ',R_TRUE

C GET THE POSITION ANGLE IN THE SKY PLANE

      SUM         = T_ANOM + ARG_PERI

      THEST       = ATAN(TAN(SUM)*COS(INCLIN))

C THEST HERE IS (OF COURSE!) REALLY THEST-ASC_NODE.  ASC_NODE 
C GETS SUBTRACTED OFF BELOW.

C     WRITE (PR,'(1X,A,F7.2,2X,F7.2,2X,F7.2)')
C    1     'TRUE2APP: ARG_PERI, SUM, THEST-ASC_NODE BEFORE CHECK: ',
C    2      ARG_PERI/RPDEG, SUM/RPDEG, THEST/RPDEG

C MAKE SURE THAT (THEST-ASC_NODE) AND (T_ANOM+ARG_PERI) END UP
C ON THE SAME BRANCH OF ARCTAN.

      IF ( (COS(SUM).GE.0. .AND. COS(THEST).LT.0.)
     1                     .OR.
     2     (COS(SUM).LT.0. .AND. COS(THEST).GE.0.) ) THEN
          THEST = THEST + PI
      END IF
      IF(THEST.GE.0.) THEN
          THEST = DMOD(THEST,TWOPI)
      ELSE
          THEST = TWOPI + DMOD(THEST,TWOPI)
      END IF

C     WRITE (PR,'(1X,A,F7.2,2X,F7.2,2X,F7.2)')
C    1     'TRUE2APP: ARG_PERI, SUM, THEST-ASC_NODE AFTER CHECK:  ',
C    2      ARG_PERI/RPDEG, SUM/RPDEG, THEST/RPDEG

      THEST       = THEST + ASC_NODE

      IF(THEST.GE.0.) THEN
          THEST = DMOD(THEST,TWOPI)
      ELSE
          THEST = TWOPI + DMOD(THEST,TWOPI)
      END IF

C GET THE SEPARATION IN THE SKY PLANE

      REST        = R_TRUE * COS(T_ANOM + ARG_PERI) / 
     1                    COS(THEST - ASC_NODE)

C CALCULATE THE RATE OF CHANGE OF SEPARATION AND POSITION ANGLE

      DV_DTIME = (TWOPI/PERIOD)*(1.+E_TRUE*COS(T_ANOM))/
     1               ((1.-E_TRUE**2)**1.5)
      DTH_DTIME = ((COS(THEST - ASC_NODE))/(COS(T_ANOM+ARG_PERI)))**2 *
     1                COS(INCLIN)*DV_DTIME
      DR_DT1 = REST*E_TRUE*SIN(T_ANOM)/(1.+E_TRUE*COS(T_ANOM))
      DR_DT2 = (R_TRUE**2/(2.*REST))*
     1                SIN(2.*(T_ANOM+ARG_PERI))*
     2                (SIN(ASC_NODE)*SIN(INCLIN)-1.)
      DR_DTIME = (DR_DT1+DR_DT2)*DV_DTIME

      DTH_DHR = DTH_DTIME/(RPDEG*24.)     ! DEGREES PER HOUR
      DR_DHR = DR_DTIME/24.               ! MAS PER HOUR

D     WRITE(PR,'(1X,A,F11.3,2(2X,F7.2),2(2X,F8.5))') 
D    1         'TRUE2APP: T, M V, DTH/DT, DR/DT:  ',
D    2         T,M_ANOM/RPDEG,T_ANOM/RPDEG,
D    3         DTH_DHR,DR_DHR

C GET THE R.A. AND DEC. SEPARATIONS

      RAEST       = REST * SIN(THEST)
      DECEST      = REST * COS(THEST)

C      WRITE(PR,'(1X,A,F7.2,A,F7.2,A,E10.3,1X,E10.3)') 
C     1            ' TRUE2APP: T_ANOM+ARG_PERI=',T_ANOM+ARG_PERI,
C     2            '     THEST-ASC_NODE=',THEST-ASC_NODE,
C     3            '     REST=',REST

C EVALUATE THE DERIVATIVES
C FIRST, SOME USEFUL EXPRESSIONS

      COSV2        = COS(T_ANOM/2.)
      COSV_AP      = COS(T_ANOM+ARG_PERI)
      SINV_AP      = SIN(T_ANOM+ARG_PERI)
      COSTH_AN     = COS(THEST-ASC_NODE)
      COSI         = COS(INCLIN)
      SINI         = SIN(INCLIN)
      SINV         = SIN(T_ANOM)
      TANV_AP      = TAN(T_ANOM+ARG_PERI)

      E_COSE       = 1. - E_TRUE * COS(E_ANOM)
      E_COSV       = 1. + E_TRUE * COS(T_ANOM)
      SEC_SEC      = (COSV2/COS(E_ANOM/2.))**2
     1                    * SQRT( (1.+E_TRUE) / (1.-E_TRUE) )
      COS2_COS2SIN2= COSV_AP**2 + SINV_AP**2 * COSI**2
      COSSINSIN2   = COSV_AP * SINV_AP * SINI**2

      R_RHO        = R_TRUE/REST
      R2_RHO       = R_TRUE**2/REST

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 THE SEPARATION R IN THE TRUE ORBIT W.R.T. SEMIMAJOR 
C AXIS, ECCENTRICITY, PERIOD, AND EPOCH:

      DRDA         = R_TRUE / A_TRUE

      DRDE         = -( 2.*A_TRUE*E_TRUE
     1                  + R_TRUE * (COS(T_ANOM)-E_TRUE*SINV*DVDE) ) / 
     2               E_COSV

      DRDP         = R_TRUE * E_TRUE * SINV * DVDP / E_COSV

      DRDT         = DRDP * DVDT / DVDP

C DERIVATIVES OF THE POSITION ANGLE THETA_EST IN THE APPARENT ORBIT W.R.T. 
C SEMIMAJOR AXIS, ARGUMENT OF PERIASTRON, ECCENTRICITY, PERIOD, EPOCH, 
C INCLINATION, AND ASCENDING NODE.

      DTHETADA     = 0.

      DTHETADAP    = COSI / ( ( 1 + (TANV_AP*COSI)**2 ) * COSV_AP**2 )

      DTHETADE     = DTHETADAP * DVDE
 
      DTHETADP     = DTHETADAP * DVDP
 
      DTHETADT     = DTHETADAP * DVDT
 
      DTHETADI     = -TANV_AP * SINI / 
     1                 ( 1. + TANV_AP**2 * COSI**2 ) 

      DTHETADAN    = 1.

C DERIVATIVES OF THE SEPARATION R_EST IN THE APPARENT ORBIT W.R.T. 
C SEMIMAJOR AXIS, ECCENTRICITY, PERIOD, EPOCH, INCLINATION, 
C ASCENDING NODE, AND ARGUMENT OF PERIASTRON.

      DRHODA       = REST / A_TRUE

      DRHODE       = R_RHO * DRDE * COS2_COS2SIN2
     1                  - R2_RHO * COSSINSIN2 * DVDE

      DRHODP       = R_RHO * DRDP * COS2_COS2SIN2
     1                  - R2_RHO * COSSINSIN2 * DVDP

      DRHODT       = R_RHO * DRDT * COS2_COS2SIN2
     1                  - R2_RHO * COSSINSIN2 * DVDT

      DRHODI       = -R2_RHO * SINV_AP**2
     1                   * SINI * COSI

      DRHODAN      = 0.

      DRHODAP      = -R2_RHO * COSV_AP * SINV_AP * SINI**2

C NOW PUT THE DERIVATIVES INTO THE DYDA ARRAY:

      DYDA(1,1)    = DRHODA
      DYDA(2,1)    = DRHODE
      DYDA(3,1)    = DRHODP
      DYDA(4,1)    = DRHODT
      DYDA(5,1)    = DRHODI
      DYDA(6,1)    = DRHODAN
      DYDA(7,1)    = DRHODAP
      
      DYDA(1,2)    = DTHETADA
      DYDA(2,2)    = DTHETADE
      DYDA(3,2)    = DTHETADP
      DYDA(4,2)    = DTHETADT
      DYDA(5,2)    = DTHETADI
      DYDA(6,2)    = DTHETADAN
      DYDA(7,2)    = DTHETADAP
  
C     WRITE(PR,'(1X,F8.3,2(1X,F7.3),2(1X,F8.4),1X,F7.3,14(1X,G11.3))')
C    1    T-2440000., E_ANOM, T_ANOM, R_TRUE, REST, THEST,
C    2    DRHODA, DRHODE, DRHODP, DRHODT, DRHODI, DRHODAN, DRHODAP,
C    3    DTHETADA, DTHETADE, DTHETADP, DTHETADT, DTHETADI, DTHETADAN,
C    4    DTHETADAP

      RETURN

      END

