NavList:

TMI, here's the USNO NOVAS ephemeris routine:

 
      SUBROUTINE DLIGHT (POS1,PE,DIFLT)
*
*     THIS SUBROUTINE RETURNS THE DIFFERENCE IN LIGHT-TIME, FOR A STAR,
*     BETWEEN THE BARYCENTER OF THE SOLAR SYSTEM AND THE OBSERVER (OR
*     THE GEOCENTER).
*
*          POS1   = POSITION VECTOR OF STAR, WITH RESPECT TO ORIGIN AT
*                   SOLAR SYSTEM BARYCENTER (IN)
*          PE     = POSITION VECTOR OF OBSERVER (OR THE GEOCENTER),
*                   WITH RESPECT TO ORIGIN AT SOLAR SYSTEM BARYCENTER,
*                   COMPONENTS IN AU (IN)
*          DIFLT  = DIFFERENCE IN LIGHT TIME, IN THE SENSE STAR TO
*                   BARYCENTER MINUS STAR TO EARTH, IN DAYS (OUT)
*
*     -OR-
*
*     THIS SUBROUTINE RETURNS THE LIGHT-TIME FROM THE OBSERVER (OR THE
*     GEOCENTER) TO A POINT ON A LIGHT RAY THAT IS CLOSEST TO A
*     SPECIFIC SOLAR SYSTEM BODY.
*
*          POS1   = POSITION VECTOR TOWARD OBSERVED OBJECT, WITH RESPECT
*                   TO ORIGIN AT OBSERVER (OR THE GEOCENTER) (IN)
*          PE     = POSITION VECTOR OF SOLAR SYSTEM BODY, WITH RESPECT
*                   TO ORIGIN AT OBSERVER (OR THE GEOCENTER), COMPONENTS
*                   IN AU (IN)
*          DIFLT  = LIGHT TIME TO POINT ON LINE DEFINED BY POS1 THAT IS
*                   CLOSEST TO SOLAR SYSTEM BODY (POSITIVE IF LIGHT
*                   PASSES BODY BEFORE HITTING OBSERVER, I.E., IF
*                   POS1 IS WITHIN 90 DEGREES OF PE)(OUT)
*
*
      DOUBLE PRECISION POS1,PE,DIFLT,C,DIS,U1,DSQRT
      DIMENSION POS1(3), PE(3), U1(3)
      SAVE
 
      DATA NTIMES / 0 /
 
      NTIMES = NTIMES + 1
      IF (NTIMES.EQ.1) THEN
*         GET C, THE SPEED OF LIGHT IN AU/DAY
          CALL ASTCON ('C(AU/DAY)',1.D0,C)
      END IF
 
*     FROM POS1, FORM UNIT VECTOR U1 IN DIRECTION OF STAR OR
*     LIGHT SOURCE
      DIS = DSQRT ( POS1(1)**2 + POS1(2)**2 + POS1(3)**2 )
      DO 20 J=1,3
   20 U1(J) = POS1(J) / DIS
 
*     LIGHT-TIME RETURNED IS THE PROJECTION OF VECTOR PE ONTO THE UNIT
*     VECTOR U1 (FORMED FROM POS1), DIVIDED BY THE SPEED OF LIGHT
      DIFLT = ( PE(1)*U1(1) + PE(2)*U1(2) + PE(3)*U1(3) ) / C
 
      RETURN
 
      END