C----------------------------------------------------------------------
C
C  R A D E C
C
C  COMPUTE RIGHT ASCENSION AND DECLINATION FROM RECTANGULAR
C  HELIOCENTRIC ECLIPTIC COORDINATES
C
C----------------------------------------------------------------------
      SUBROUTINE RADEC(
     .             X,Y,Z, X0,Y0,Z0,
     .             OLDRA,OLDDEC,DAYS,J,T,
     .             RA,DEC,DIST)
      INTEGER J  
      REAL
     .    X,Y,Z ,
     .    X0,Y0,Z0, 
     .    OLDRA,OLDDEC, 
     .    DAYS, 
     .    T,    
     .    RA,   
     .    DEC,  
     .    DIST  
C
C  CONSTANTS
C
      REAL  PI,PI2,CPI,DTORAD,RTODEG,
     .                  GRAVIT,OBLIQ
      COMMON /CONST/    PI,PI2,CPI,DTORAD,RTODEG,
     .                  GRAVIT,OBLIQ
C 
C  LOCAL VARIABLES
C
      REAL
     .    U,K,D,G,  
     .    DRA,DDEC,TAU,
     .    DPSI,DEPS,
     .    XL,YL,ZL  
C
      XL=X-X0
      YL=Y-Y0
      ZL=Z-Z0
      DIST=SQRT(XL*XL+YL*YL+ZL*ZL)
      DEC=ASIN(ZL/DIST)
      RA=ATAN2(YL,XL)
      IF (RA .LT. 0.0) RA=RA+PI2
C
C  ANNUAL  ABERRATION CORRECTION
C
      IF (J .NE. 1) THEN
         D=ZL/DIST
         G=SQRT(1-D*D)
         U=-X0/SQRT(X0*X0+Y0*Y0+Z0*Z0)
         K=SQRT(1.0E+0-U*U)
         DRA=1.363*(U*COS(OBLIQ)*COS(RA)+K*SIN(RA))/G
         DDEC=20.445*(U*(G*SIN(OBLIQ)-D*COS(OBLIQ)*SIN(RA))
     .        +D*K*COS(RA))
         RA=RA-DRA*15/3600*DTORAD
         DEC=DEC-DDEC/3600*DTORAD
      ENDIF
C
C  PLANETARY ABERRATION 
C
      D=DIST/173.1446577 
      OLDRA=OLDRA*DTORAD*15
      OLDDEC=OLDDEC*DTORAD
      IF ((OLDRA .NE. 0.0) .OR. (OLDDEC .NE. 0.0)) THEN
         DRA=RA-OLDRA
         IF (DRA/DTORAD/15 .LT. -12.0)  DRA=DRA+PI2
         IF (DRA/DTORAD/15 .GT. +12.0)  DRA=DRA-PI2
         DRA=(DRA)*D/DAYS
         DDEC=(DEC-OLDDEC)*D/DAYS
         RA=RA-DRA
         DEC=DEC-DDEC
      ENDIF
C
C  NUTATION
C  AE S.B20 (1982)
C
      D=T+(2415020.0E+0-2444969.5E+0)
      DPSI=(-0.0048*SIN((113.2-0.053*D)*DTORAD)
     .          -0.0004*SIN((198.7+1.971*D)*DTORAD))*DTORAD
      DEPS=( 0.0026*COS((113.2-0.053*D)*DTORAD)
     .         +0.0002*COS((198.7+1.971*D)*DTORAD))*DTORAD
      DRA=(COS(OBLIQ)+SIN(OBLIQ)*SIN(RA)*TAN(DEC))*DPSI
     .     -COS(RA)*TAN(DEC)*DEPS
      DDEC=SIN(OBLIQ)*COS(RA)*DPSI+SIN(RA)*DEPS
      RA=RA+DRA
      DEC=DEC+DDEC
C
      RA=RA*RTODEG/15   
      DEC=DEC*RTODEG    
      RETURN
      END
