
>>>>>>>>>>>>>>>>>>> FILE: ../comp/cadget.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE CADINI( TREF, NC, JGS, SBP, BPL)
C
      USE DERIVS_MOD, ONLY: DERIVS_CA
      USE PROB_IP_MOD, ONLY: PROB_IP, PROB_IP_MSG
C
      DOUBLE PRECISION   TREF
      INTEGER            NC
      LOGICAL            JGS, SBP, BPL
C
C     SUBROUTINE CADGET( RI, NP, ICAB, CAD)
      DOUBLE PRECISION   RI(*), CAD(11)
      INTEGER            NP, ICAB
C
C     SUBROUTINE CADGAD( CAA, P, PX)
      DOUBLE PRECISION   CAA(9), P(*), PX(6,*)
C
C     SUBROUTINE CADSTA( XP, MU, RP)
      DOUBLE PRECISION   XP(6), MU, RP
C
C     SUBROUTINE CADPRT( ICAB, CAD)
C
C     SUBROUTINE CADWRT( U, ICAB, CAD, CAA, P, PX)
      INTEGER            U
C
C-----------------------------------------------------------------------
C  CADINI (Close Approach Data INItialize) initializes this package.
C
C  Inputs:
C   TREF   Reference time (JED).
C   NC     Number of columns being integrated: NC = 1 if partials are
C          not being integrated; NC = NDP + 1 if partials are being
C          integrated, where NDP is the number of dynamic parameters.
C   JGS    Jupiter/Galilean Satellites separation flag:
C          TRUE to separate them, FALSE to combine them.
C   SBP    Small bodies perturbations flag: TRUE to include them.
C   BPL    B-Plane flag: TRUE to use b-plane as target plane.
C
C-----------------------------------------------------------------------
C  CADGET (Close Approach Data GET) computes close approach data for the
C  next close approach in the current integration step.  CADGET should
C  be called once for each close approach occurring in the current step.
C  The routine uses one of two methods to compute the close approach
C  data: if BPL is TRUE in the call to CADINI, the b-plane is used as
C  the target plane; otherwise, the plane perpendicular to the relative
C  velocity vector at close approach is used as the target plane.  If
C  the relative motion is not hyperbolic, the latter plane is used
C  regardless of the setting of BPL.
C
C  Inputs:
C   RI(1..NP*(NP+1)/2)  The epoch sqrt covariance matrix, upper-triangu-
C                       lar, vector-stored.
C   NP                  Number of parameters.
C
C  Outputs:
C   ICAB        Index of Close Approach Body.
C   CAD(1..11)  Close approach data, as follows:
C                CAD(1):  Time of CA (days relative to ref time).
C                CAD(2):  Nominal close approach distance (AU).
C                CAD(3):  Relative velocity at close approach (km/s).
C                CAD(4):  Minimum 3-sigma close approach distance (AU).
C                CAD(5):  Maximum 3-sigma close approach distance (AU).
C                CAD(6):  Uncertainty in CA time, 3-sigma (min).
C                CAD(7):  Target plane ellipse semi-major axis (km).
C                CAD(8):  Target plane ellipse semi-minor axis (km).
C                CAD(9):  Target plane ellipse angle from range (deg).
C                CAD(10): Minimum number of sigmas for impact.
C                CAD(11): Probability of impact.
C
C-----------------------------------------------------------------------
C  CADGAD (Close Approach Get Additional Data) returns additional data
C  associated with the close approach.  The close approach state is
C  represented either in scaled b-plane elements or target-plane-frame
C  cartesian position and velocity.
C
C  Outputs:
C   CAA(1..9)   Close approach additional data, as follows:
C                CAA(1):    Position angle of center of ellipse, in
C                           target plane, measured counterclockwise from
C                           projection of inertial +z-axis (deg).
C                CAA(2..4): Inertial-frame components of unit vector
C                           into target plane.
C                CAA(5):    Scaled b magnitude, if b-plane state (km).
C                CAA(6):    T_lin - T_CA (s).
C                CAA(7):    V_infinity, if b-plane state (km/s);
C                           zero otherwise.
C                CAA(8..9): Change in target plane coordinates to get
C                           impact with minimum number of sigmas (km).
C   P(1..NP*(NP+1)/2)  Close approach state covariance matrix, upper-
C                      triangular, vector-stored.  The first 6 rows and
C                      columns are for close approach state; the
C                      remaining are for the other estimated parameters.
C   PX(1..6,1..NDP)    Partials of close approach state with respect to
C                      epoch state and the other estimated parameters.
C
C-----------------------------------------------------------------------
C  CADSTA (Close Approach STAte) returns the planetocentric state at 
C  closest approach. It also returns the planets size and mass.
C
C  Outputs:
C   XP(1..6 )   Planetocentric close approach state
C   MU          GM of planet
C   RP          Radius of planet
C
C-----------------------------------------------------------------------
C  CADPRT (Close Approach Data PRinT) prints a row of the close approach
C  table on unit 6.  If this is the first call to CADPRT, a header is
C  printed.  If the b-plane method was not used to compute the data,
C  an asterisk flag is prepended to the body name in the table.
C
C  Inputs:
C   CAD(1..11)  Close approach data (same as above).
C
C-----------------------------------------------------------------------
C  CADWRT (Close Approach Data WRiTe) writes the close approach data to
C  in formatted form to the specified unit number.
C
C  Inputs:
C   U                  Unit number.
C   ICAB               Index of Close Approach Body.
C   CAD(1..11)         Close approach data (same as above).
C   CAA(1..9)          Close approach additional data (same as above).
C   P(1..NP*(NP+1)/2)  Close approach state covariance matrix, upper-
C                      triangular, vector-stored (same as above).
C   PX(1..6,1..NDP)    Partials of close approach state with respect to
C                      epoch state and the other estimated parameters.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/etmut.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE  ETMUT( JD, DT, ERR)
      DOUBLE PRECISION   JD, DT
      INTEGER            ERR
C
C     SUBROUTINE ETMINI( U, DEMN)
      INTEGER            U
      REAL               DEMN
C
C-----------------------------------------------------------------------
C  ETMUT (ET Minus UT) computes the value of ET-UTC.  This version is
C  based on L.V.Morrison's values of ET-UTC.  ETMINI is used to
C  initialize the post-1972 table.
C
C  Input:
C   JD     UTC Time (Julian Date).
C
C  Output:
C   DT     ET-UTC (s).
C   ERR    Error indicator:  0 = no error
C                            1 = beyond table data
C                            2 = before year -3002
C
C-----------------------------------------------------------------------
C  ETMINI (ET Minus ut INItialize) opens the ET-UT data file, reads and
C  initializes the post-1972 table, and closes the file.
C  For portability reasons, the user must supply a routine which
C  performs the OPEN on the file, and sets the I/O unit number; the
C  calling sequence to this routine must be as follows:
C      SUBROUTINE OPNETM( UNIT )
C      INTEGER            UNIT
C  where UNIT is the I/O unit number, an input of the routine.
C
C  Inputs:
C   U     The I/O unit number to be passed to OPNETM.
C   DEMN  Derivative of Moon's mean motion, ie, lunar tidal accelera-
C         tion (arc-sec/century^2).  Some possible values:
C           -26.        Morrison, MNRAS 187,41 (1979)
C           -22.44      Spencer Jones (Clemence), MNRAS 99,541 (1939)
C                        (GMC) AJ 53,169 (1948)
C           -23.8 +- 4  Williams, Sinclair, Yoder,
C                        Geophys.Rev.Lett. 5, 943 (1978)
C
C-----------------------------------------------------------------------
C Notes:
C  VALUES 1621 TO 1860 ARE PRIV. COMM. MORRISON (1980)
C  VALUES 1861 TO 1976 ARE MORRISON GJRAS 58, 349 (1979)
C  VALUES POST 1976 ARE BIH CIRC. D + 32.184
C       AND ARE SMOOTHED VIA MORRISON FORMULA
C         35 S(0)= 17 F(0) + 12( F(1)+F(-1) ) -3 ( F(2)+F(-2) )
C
C  LUNAR LONGITUDE IS L= L0 + L1*T + DEMON/2 *T**2
C-- NOTE THAT DEMON IS NDOT AND TWICE THE COEFFICIENT OF T**2
C   IN LONGITUDE OF MOON.
C-- IF ALTERNATE VALUE OF DEMON IS USED, THEN THE ADDITION TO
C   MORRISON VALUES (PRIOR TO 1955.5 ONLY) IS
C
C           D-DELTA T (NEW - MORRISON) = -0.911 * (DEMON+26) *T**2
C                 WHERE D-DELTA T IS IN SEC OF TIME AND
C                 DEMON IS ARCSEC PER CENTURY-SQUARED WITH
C                 T = (YEAR - 1955.5) / 100
C
C-- I.E., IF USE S-J VALUE NDOT=-22.44, GET D-DELTA T = -3.27 *T**2 SEC
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/nutang.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE NUTANG( TD, NUT)
      DOUBLE PRECISION   TD, NUT(4)
C
C-----------------------------------------------------------------------
C  NUTANG (NUTation ANGles) computes the two nutation angles (in
C  longitude and obliquity) and their rates, for a specified Julian
C  date, using the 1980 IAU Theory of Nutation.  Developed by J. Wahr,
C  this theory is comprised of two 106-term series.
C
C  Input:
C   TD         Julian date (JED).
C
C  Output:
C   NUT(1..4)  Nutation angles and their rates:
C               NUT(1) is the nutation in longitude (rad).
C               NUT(2) is the nutation in obliquity (rad).
C               NUT(3) is the nutation rate in longitude (rad/d).
C               NUT(4) is the nutation rate in obliquity (rad/d).
C
C  References:
C    1. Supplement to the 1984 Astronomical Almanac, pp. S21-S26.
C    2. P.K. Seidelmann et al., "1980 Theory of Nutation",
C          Celest. Mech., 27, pp. 79-105, 1982.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/nutatd.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE NUTATD( OBM, NUT, RM, RT)
C
C     SUBROUTINE NUTAMD( OBM, NUT, RT, RM)
C
C     SUBROUTINE NUTVMD( OBMD, NUTD, RM, VT, VM)
C
      DOUBLE PRECISION   OBM, NUT(2), RM(3), RT(3)
      DOUBLE PRECISION   OBMD, NUTD(2), VT(3), VM(3)
C
C-----------------------------------------------------------------------
C  NUTATD (NUTAte to True of Date) converts the components of a position
C  vector from the system with mean equator and equinox of date to the
C  system with true equator and equinox of date.  This transformation
C  represents the effects of nutation.
C
C  Inputs:
C   OBM       Mean Obliquity of the ecliptic at date (rad).
C   NUT(1..2) Nutation angles at date:
C              NUT(1) is the nutation in longitude (rad),
C              NUT(2) is the nutation in obliquity (rad).
C   RM(1..3)  Position vector components in mean-of-date system.
C
C  Output:
C   RT(1..3)  Position vector components in true-of-date system.
C             RT can overwrite RM in the calling routine.
C
C-----------------------------------------------------------------------
C  NUTAMD (NUTate to Mean of Date) converts the components of a position
C  vector from the system with true equator and equinox of date to the
C  system with mean equator and equinox of date.  This transformation
C  represents the effects of nutation, and is the inverse of the trans-
C  formation computed by entry point NUTATD.  The arguments are
C  analogous:
C
C  Inputs:
C   OBM       Mean Obliquity of the ecliptic at date (rad).
C   NUT(1..2) Nutation angles at date:
C              NUT(1) is the nutation in longitude (rad),
C              NUT(2) is the nutation in obliquity (rad).
C   RT(1..3)  Position vector components in true-of-date system.
C
C  Output:
C   RM(1..3)  Position vector components in mean-of-date system.
C             RM can overwrite RT in the calling routine.
C
C-----------------------------------------------------------------------
C  NUTVMD (NUTate Velocity to Mean of Date) converts the components of
C  a velocity vector from the true-of-date system to the mean-of-date
C  system.  This entry point is an extension of NUTAMD, which must be
C  called first. The date TD is assumed to be the same as that of the
C  previous call to NUTAMD, and the output RM from that previous call to
C  NUTAMD is an input to this entry point.
C
C  Inputs:
C   OBMD       Rate of change of mean obliquity (rad/d).
C   NUTD(1..2) Nutation angle rates at time TD:
C               NUTD(1) is the nutation rate in longitude (rad/d),
C               NUTD(2) is the nutation rate in obliquity (rad/d).
C   RM(1..3)   Position vector components in mean-of-date system.
C   VT(1..3)   Velocity vector components in true-of-date system.
C
C  Output:
C   VM(1..3)   Velocity vector components in mean-of-date system.
C              VM can overwrite VT in the calling routine.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/objsta.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE OBJSTA( FWD, TO, TI, XI, RS, TC, X, RHO, SUCC)
      LOGICAL FWD
      DOUBLE PRECISION   TO, TI, XI(3), RS(3), TC, X(6), RHO
      LOGICAL SUCC
C
C     SUBROUTINE OBJSDT( DELT)
      DOUBLE PRECISION   DELT
C
C-----------------------------------------------------------------------
C  OBJSTA (OBJect STAte) computes the state (inertial-frame position and
C  velocity components) of the object at the time the light left it in
C  order to be observed at the station at time TO.  The state is
C  obtained via interpolation within the current integration step.  If
C  the time the light left the object is not within the current step,
C  the routine returns with SUCC=.FALSE. and does not compute the state.
C
C  Inputs:
C   FWD       Logical to indicate whether the integration is moving 
C             forward or backward.
C   TO        Time of observation (ET days wrt reference time).
C   TI        Current integration time (ET days wrt reference time).
C   XI(1..3)  Inertial-frame position components of object wrt Sun at
C             time TI (AU).
C   RS(1..3)  Inertial-frame position components of station wrt Sun at
C             time TO (AU).
C
C  Outputs:
C   TC        Time the light left the object (ET days wrt reference
C             time).
C   X(1..6)   Inertial-frame position and velocity components of object
C             wrt Sun at time TC (AU, AU/d); not computed if SUCC=.FALSE.
C   RHO       Distance travelled by the light (AU); not computed if
C             SUCC=.FALSE.
C   SUCC      Logical to indicate success.
C
C-----------------------------------------------------------------------
C  OBJSDT (OBJect State Delta T) initializes the value of DT. If not 
C  initialized then the previously computed value will be used (or zero 
C  on the first call).
C
C  Input:
C   DELT      Estimated one-way light time from object to observer.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/obliqm.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE OBLIQM( TD, OBM)
      DOUBLE PRECISION   TD, OBM(2)
C
C     SUBROUTINE OBLSYS( J2000, TREF)
      LOGICAL            J2000
      DOUBLE PRECISION   TREF
C
C-----------------------------------------------------------------------
C  OBLIQM (OBLIQuity, Mean) computes the mean obiquity of the ecliptic
C  at a specified date TD, and its rate of change.  The mean ecliptic
C  can be either that defined by the FK4/B1950 system of astronomical
C  constants, or the FK5/J2000 system of constants, as specified in the
C  call to entry point OBLSYS.  If OBLSYS is not called, the system
C  defaults to FK4/B1950.  The date TD can be either a full Julian date,
C  or a time interval in days from a reference time TREF, which is also
C  specified in the initialization call to OBLSYS.  If TD is a full
C  Julian date, TREF should be 0.D0 in the call to OBLSYS.
C
C  Input:
C   TD         Date of the true-of-date system, in days relative to the
C              reference time TREF specified in the call to NUTSYS.
C
C  Output:
C   OBM(1..2)  Mean obliquity of the ecliptic and its rate at time TD:
C               OBM(1) is the mean obliquity (rad).
C               OBM(2) is the rate of change of mean obliquity (rad/d).
C
C
C  References:
C    1. Explanatory Supplement to the Astronomical Ephemeris,
C         H.M. Nautical Almanac Office, 1961.
C    2. Lieske, J., "Expressions for the Precession Quantities and
C         Their Partial Derivatives", JPL Tech Rep. 32-1044, 1967.
C    3. Lieske, J., "Precession Matrix Based on IAU (1976) System of
C         Astronomical Constants", Astron. Astrophys. 73, 282-284, 1979.
C    4. The Astronomical Almanac, 1984, U.S. Government Printing Office.
C
C-----------------------------------------------------------------------
C  OBLSYS (OBLiquity SYStem) defines the system of astronomical cons-
C  tants to be used in computing the obliquity (FK4/B1950 or FK5/J2000),
C  and sets TREF, which is an arbitrary reference time from which TD is
C  measured.  If TD will be a full Julian date, TREF should be 0.D0.
C  OBLSYS should be called before OBLIQM is called.  If OBLSYS is not
C  called, the system defaults to FK4/B1950 and TREF defaults to 0.D0.
C
C  Inputs:
C   J2000  Set TRUE to indicate the FK5/J2000 system is to be used;
C          set FALSE to indicate the FK4/B1950 system.
C   TREF   Reference time for TD (JED).
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/precmd.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PRECMD( TD, RI, RM)
C
C     SUBROUTINE PRECIN( TD, RM, RI)
C
C     SUBROUTINE PREVIN( RI, VM, VI)
C
      DOUBLE PRECISION   TD, RI(3), RM(3), VM(3), VI(3)
C
C     SUBROUTINE PRESYS( J2000, TREF)
      LOGICAL            J2000
      DOUBLE PRECISION   TREF
C
C-----------------------------------------------------------------------
C  PRECMD (PRECess to Mean of Date) converts the components of a
C  position vector from the inertial reference system to the system with
C  mean equator and equinox of date.  This transformation represents the
C  effects of precession.  The input reference system can be either
C  FK4/B1950 or FK5/J2000, as specified in an initialization call to
C  routine PRESYS.  The date TD can be either a full Julian date, or a
C  time interval in days from a reference time TREF, which is also
C  specified in the initialization call to PRESYS.  If TD is a full
C  Julian date, TREF should be 0.D0 in the call to PRESYS.
C
C  Inputs:
C   TD        Date of the mean-of-date system, in days relative to the
C             reference time TREF specified in the call to PRESYS.
C   RI(1..3)  Position vector components in inertial system.
C
C  Output:
C   RM(1..3)  Position vector components in mean-of-date system.
C             RM can overwrite RI in the calling routine.
C
C  References:
C    1. Explanatory Supplement to the Astronomical Ephemeris,
C         H.M. Nautical Almanac Office, 1961.
C    2. Lieske, J., "Expressions for the Precession Quantities and
C         Their Partial Derivatives", JPL Tech Rep. 32-1044, 1967.
C    3. Lieske, J., "Precession Matrix Based on IAU (1976) System of
C         Astronomical Constants", Astron. Astrophys. 73, 282-284, 1979.
C
C-----------------------------------------------------------------------
C  PRECIN (PRECess to INertial) converts the components of a position
C  vector from the system with mean equator and equinox of date to the
C  inertial system.  This transformation represents the effects of
C  precession, and is the inverse of the transformation computed by
C  PRECMD.  The output system can be either FK4/B1950 or FK5/J2000, as
C  specified in an initialization call to PRESYS.  The arguments are
C  analogous to those of PRECMD:
C
C  Inputs:
C   TD        Date of the mean-of-date system, in days relative to the
C             reference time TREF specified in the call to PRESYS.
C   RM(1..3)  Position vector components in mean-of-date system.
C
C  Output:
C   RI(1..3)  Position vector components in inertial system.
C             RI can overwrite RM in the calling routine.
C
C-----------------------------------------------------------------------
C  PREVIN (PREcess Velocity to INertial) converts the components of a
C  velocity vector from the mean-of-date system to the inertial system.
C  This entry point is an extension of PRECIN, which must be called
C  first.  The date TD is assumed to be the same as that of the previous
C  call to PRECIN, and the output RI from that previous call to PRECIN
C  is an input to this entry point.
C
C  Inputs:
C   RI(1..3)  Position vector components in inertial system.
C   VM(1..3)  Velocity vector components in mean-of-date system.
C
C  Output:
C   VI(1..3)  Velocity vector components in inertial system.
C             VI can overwrite VM in the calling routine.
C
C-----------------------------------------------------------------------
C  PRESYS (PREces SYStem) sets the system (FK4/B1950 or FK5/J2000) used
C  by PRECMD, and sets TREF, which is an arbitrary reference time from
C  which TD is measured.  If TD will be a full Julian date, TREF should
C  be 0.D0.  PRESYS should be called before PRECMD is called.  If PRESYS
C  is not called, the system defaults to FK4/B1950 and TREF defaults to
C  0.D0.
C
C  Inputs:
C   J2000  Set TRUE to indicate the FK5/J2000 system is to be used;
C          set FALSE to indicate the FK4/B1950 system.
C   TREF   Reference time for TD (JED).
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/radec.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE RADEC ( R, RT, TO, RQPAR, RA, DEC, PA, PD)
      DOUBLE PRECISION   R(3), RT(3), TO, RA, DEC, PA(4), PD(4)
      LOGICAL            RQPAR
C
C     SUBROUTINE RADCLM( CLS0, CLN, TCL, ESTCL)
      DOUBLE PRECISION   CLS0, CLN, TCL
      LOGICAL            ESTCL
C
C     SUBROUTINE RADCLU( DS0, CLS0)
      DOUBLE PRECISION   DS0
C
C     SUBROUTINE RADGCP( CLS0, CLN, TCL)
C
C-----------------------------------------------------------------------
C  RADEC (Right Ascension and DEClination) computes the right ascension
C  and declination of an object, given its heliocentric position at the
C  time the light left the object, and the heliocentric position of the
C  observer at the time the object was observed.  If requested, the
C  routine applies a center-of-light offset model in which the observed
C  center of light is offset from the center of mass radially towards
C  the Sun by a distance which varies according to an inverse power of
C  the distance from the Sun.  If requested, the routine also computes
C  the partial derivatives of the right ascension and declination with
C  respect to the position of the object at the time the light left it.
C  The partials are requested via a flag in the calling sequence.  The
C  center of light offset model is set up by an earlier call to RADCLM.
C
C  Inputs:
C   R(1..3)   Inertial-frame heliocentric position components of the
C             object at the time the light left it (AU).
C   RT(1..3)  Inertial-frame heliocentric position components of the
C             observer at the time of the observation (AU).
C   TO        Time of the observation (d wrt ref time).
C   RQPAR     Set TRUE to request computation of partials.
C
C  Outputs:
C   RA        Computed right ascension (rad); range -pi..pi.
C   DEC       Computed declination (rad); range -pi/2..pi/2.
C   PA(1..4)  Partials of right ascension w.r.t. the position components
C             R, output only if RQPAR is TRUE.  If the center-of-light
C             model is enabled, PA(4) is the partial w.r.t. S0.
C   PD(1..4)  Partials of declination w.r.t. the position components R,
C             output only if RQPAR is TRUE.  If the center-of-light
C             model is enabled, PA(4) is the partial w.r.t. S0.
C
C  Note:  The computed R.A. and Dec are relative to the inertial frame
C         used to express the inputs R and RT.
C
C-----------------------------------------------------------------------
C  RADCLM (R.A. & Dec. Center of Light Model) sets up the center of
C  light model, which is otherwise disabled.  The model is enabled if
C  either S0 is nonzero or ESTCL is true.
C
C  Inputs:
C   CLS0   Center-of-light offset at 1 AU, S0 (km).
C   CLN    Inverse power by which the center-of-light offset varies.
C          Eg, CLN=2 for inverse square variation.
C   TCL    Start time for center-of-light offset model (d wrt ref time).
C   ESTCL  Set TRUE if center-of-light offset is being estimated.
C
C-----------------------------------------------------------------------
C  RADCLU (R.A. & Dec. Center of Light model Update) updates the center-
C  of-light offset by adding a correction to it.
C
C  Input:
C   DS0    Correction to S0 (km).
C
C  Output:
C   CLS0   Corrected value of S0 (km).
C
C-----------------------------------------------------------------------
C  RADGCP (R.A. & Dec. Get Center of light model Parameters) returns the
C  value of the center-of-light offset.
C
C  Outputs:
C   CLS0   Center-of-light offset at 1 AU, S0 (km).  A value of 0.D0
C           indicates that the c/l offset model is not active, in which
C           case the next two outputs are not set.
C   CLN    Inverse power by which the center-of-light offset varies.
C   TCL    Start time for center-of-light offset model (d wrt ref time).
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/sdtime.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE SDTIME( TUT1, GMST)
      DOUBLE PRECISION   TUT1, GMST
C
C     SUBROUTINE SDTIMD( OME)
      DOUBLE PRECISION   OME
C
C     SUBROUTINE SDTSYS( J2000, TREF)
      LOGICAL            J2000
      DOUBLE PRECISION   TREF
C
C-----------------------------------------------------------------------
C  SDTIME (SiDereal TIME) computes the Greenwich mean sidereal time, as
C  an angle, at a specified UT1 time.  The expression used to model the
C  rotation of the Earth is either that of the FK4/B1950 system of
C  astrodynamical constants, or the FK5/J2000 system, as specified in
C  an initialization call to entry point SDTSYS.  The input UT1 time
C  can be either a full Julian date, or a time offset in days from a
C  reference time TREF, which is also specified in the call to SDTSYS.
C  If the input time is an absolute Julian date, TREF should be
C  specified as 0.D0 in the call to SDTSYS.
C
C  Input:
C   TUT1   Universal Time (UT1), in days relative to ref time TREF.
C
C  Output:
C   GMST   Greenwich Mean Sidereal Time (rad).
C
C  References:
C    1. Explanatory Supplement to the Astronomical Ephemeris,
C         H.M. Nautical Almanac Office, 1961.
C    2. The Astronomical Almanac, 1984, U.S. Government Printing Office.
C
C-----------------------------------------------------------------------
C  SDTIMD (SiDereal TIMe Dot) computes the rate of change of sidereal
C  time at the time used in the most recent call to SDTIME.
C
C  Output:
C   OME    Angular rate of the Earth (rad/d).
C
C-----------------------------------------------------------------------
C  SDTSYS (SiDereal Time SYStem) sets the system of astrodynamical
C  constants (FK4/B1950 or FK5/J2000) to be used in the computation of
C  sidereal time, and sets TREF, which is an arbitrary reference time
C  from which TUT1 is measured.  If TUT1 will be a full Julian date,
C  TREF should be specified as 0.D0.  SDTSYS should be called before
C  SDTIME is called: if SDTSYS is not called, the system defaults to
C  FK4/B1950, and TREF defaults to 0.D0.
C
C  Input:
C   J2000     Set TRUE  to use the FK5/J2000 system;
C             set FALSE to use the FK4/B1950 system.
C   TREF      Reference time for TUT1 (JD, UT1).
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../comp/eop_mod.f90 <<<<<<<<<<<<<<<<<<<

module eop_mod

use constants_mod
use get_unit_mod

implicit none
save
private

! Public procedures
public :: eop_ini, eop_get, eop_def

! Public variables
public :: eop_msg, eoplbl, eopfng, eoput1, eoptyp, eoptim, eoptrf, eopcrf

! Public module variables
character (len=256) :: eop_msg     ! Diagnostic message
character (len=80)  :: eoplbl= " " ! EOP Label
character (len=80)  :: eopfng= " " ! EOP Fingerprint
character (len=6)   :: eoput1= " " ! EOP UT1 type (UT1 or UT1R)
character (len=6)   :: eoptyp= " " ! EOP Type
character (len=25)  :: eoptim= " " ! EOP Time of file creation
character (len=6)   :: eoptrf= " " ! EOP Terrestrial Ref. Frame
character (len=6)   :: eopcrf= " " ! EOP Celestial Ref. Frame

!===============================================================================
! The eop_mod module is used to load and interpolate Earth Orientation 
! Parameter (EOP) data from a JPL-style EOP file.
!
! Public routines: 
!    eop_ini        - Loads the EOP data into memory
!    eop_get        - Interpolates the EOP data
!    eop_def        - VLDEFs the EOPFILE VARLIST variable
!
! Public variables:
!    eop_msg     - Diagnostic message in case of error.
!    eoplbl      - EOP Label
!    eopfng      - EOP Fingerprint
!    eoput1      - EOP UT1 type (UT1 or UT1R)
!    eoptyp      - EOP Type
!    eoptim      - EOP Time of file creation
!    eoptrf      - EOP Terrestrial Ref. Frame
!    eopcrf      - EOP Celestial Ref. Frame
!
! References: 
!     http://eis.jpl.nasa.gov/nav/eop/eop.html
!     JPL-IOM 335.1-11-93  May 21, 1993 "Earth Orientation Parameter (EOP) 
!         file description and usage," Folkner, Steppe and Oliveau
!     JPL-IOM 335.5-93-666 Apr 29, 1993 "Interpolation Routine UTPM," Newhall
!     JPL-EM  314-558      Feb 26, 1993 "Frame Tie Rotation and Nutation 
!         Corrections for the ODP," Moyer
!-------------------------------------------------------------------------------
! Written       Jan 31, 2002       Steve Chesley
! Modified      May 28, 2002       Steve Chesley
!    - Handle blank lines in EOP file.
! Modified      Jan 06, 2005       Steve Chesley
!    - Fix bug related to comment skipping in eop_ini.
!
! $Id: eop_mod.f90,v 1.7 2005/01/06 20:47:07 chesley Exp $
!$==============================================================================

!^########################## EOP_INI ###########################################
subroutine eop_ini(t_ref, errcod)

real (kind=wp), intent(in) :: t_ref
integer, intent(out) :: errcod

!===============================================================================
! eop_ini loads the EOP file from the first found among:
!         1) EOPFILE variable from VARLIST
!         2) EOPFILE variable from environment
!         3) $CAET_DATA/eop.dat
! The routine also makes available the following labels and flags through 
! public character variables:
!         EOPLBL - EOP Label
!         EOPFNG - EOP Fingerprint
!         EOPUT1 - EOP UT1 type (UT1 or UT1R)
!         EOPTYP - EOP Type
!         EOPTIM - EOP Time of file creation
!         EOPTRF - EOP Terrestrial Ref. Frame
!         EOPCRF - EOP Celestial Ref. Frame
!
! Inputs:
!    t_ref  - Reference time, days relative to JD=0.
!
! Outputs:
!    errcod - Zero indicates success.
!$==============================================================================

!^########################## EOP_GET ###########################################
subroutine eop_get(utc, rates, taimutc, taimut1, p_xy, delta_nut, errcod)

real (kind=wp), intent(in) :: utc
logical, intent(in) :: rates
real (kind=wp), intent(out) :: taimutc
real (kind=wp), intent(out) :: taimut1
real (kind=wp), intent(out) :: p_xy(:)
real (kind=wp), intent(out) :: delta_nut(:)
integer, intent(out) :: errcod

!===============================================================================
! eop_get reads the internal EOP array and interpolates to the requested
! UTC time. It uses a slightly altered version of the original UTPM routine
! by Skip Newhall. The nutation rate corrections are optionally computed.
!
! Inputs:
!    utc       - UTC in days relative to t_ref for which the earth orientation
!                parameters are requested.
!    rates     - flag to determine whether nutation rates should be computed.
!
! Outputs:
!    taimutc   - TAI - UTC (sec).
!    taimut1   - TAI - UT1 (sec). [NOTE: This is _not_ TAI - UT1R.]
!    p_xy      - Polar motion X and Y calibrations (arcsec).
!    delta_nut - Nutation corrections. First two elements are [d_psi, d_eps].
!                If rates are requested then next two elements are the
!                corresponding time derivatives.
!    errcod - Zero indicates success.
!$==============================================================================

!^########################## EOP_DEF ###########################################
subroutine eop_def

!===============================================================================
! VLDEF the EOPFILE variable.
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../comp/et_minus_tai_mod.f90 <<<<<<<<<<<<<<<<<<<

module et_minus_tai_mod

use constants_mod

implicit none
save
private

! Public procedures
public :: et_minus_tai, et_minus_tai_ini

! Public variables
public :: et_minus_tai_msg

! Public module variables
character (len=256) :: et_minus_tai_msg ! Diagnostic message

!===============================================================================
! The et_minus_tai_mod module is used to compute ET - TAI.
!
! Public routines: 
!    et_minus_tai_ini - Initializes a few module parameters
!    et_minus_tai     - Returns ET-TAI at a given location
!
! Public variables:
!    et_minus_tai_msg - Diagnostic message in case of error.
!
! References: Theodore D. Moyer, "Formulation for Observed and Computed 
!    Values of Deep Space Network Data Types for Navigation" 
!    JPL Publication 00-7 (2000).
!-------------------------------------------------------------------------------
! Written       Feb 04, 2002       Steve Chesley
! Modified      Feb 20, 2002       Steve Chesley
!   - Fixed bug relating to v_light units.
! Modified      Sep 17, 2003       Steve Chesley
!   - Tweaked documentation
!
! $Id: et_minus_tai_mod.f90,v 1.4 2004/12/02 00:37:15 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../comp/stn_pos_mod.f90 <<<<<<<<<<<<<<<<<<<

module stn_pos_mod

use constants_mod
use get_unit_mod
use eop_mod
use et_minus_tai_mod

implicit none
save
private

! Public procedures
public :: stn_pos_ini, stn_pos

! Public variables
public :: stn_pos_msg

! Public module variables
character (len=256) :: stn_pos_msg ! Diagnostic message

!===============================================================================
! The stn_pos_mod module is used to compute earth-fixed station positions 
! in the planetary ephemeris frame.
!
! Public routines: 
!    stn_pos_ini - Initializes the reference time and reference frame routines.
!    stn_pos     - Computes the station position.
!
! Public variables:
!    stn_pos_msg - Diagnostic message in case of error.
!
!-------------------------------------------------------------------------------
! Based on RDRSTN/DELTIM, written by Paul Chodas in 1992.
!
! Written       Feb 01, 2002       Steve Chesley
!
! $Id: stn_pos_mod.f90,v 1.4 2004/12/02 00:37:15 chesley Exp $
!$==============================================================================

!^########################## STN_POS ###########################################
subroutine stn_pos(utc, stn_coords, radar, use_eop, compute_et, et, etmutc, &
   x_stn, x_earth, r_stn_sun, errcod)

real (kind=wp), intent(in)    :: utc
real (kind=wp), intent(in)    :: stn_coords(3)
logical,        intent(in)    :: radar
logical,        intent(in)    :: use_eop
logical,        intent(in)    :: compute_et
real (kind=wp), intent(inout) :: et
real (kind=wp), intent(out)   :: etmutc
real (kind=wp), intent(out)   :: x_stn(6)
real (kind=wp), intent(out)   :: x_earth(6)
real (kind=wp), intent(out)   :: r_stn_sun(3)
integer,        intent(out)   :: errcod

!===============================================================================
! stn_pos computes the position of an observing station with respect to
! both the Sun and the Solar System barycenter at a specified UTC time.
! Precession, nutation, polar motion and the variability of the Earth
! rotation are taken into account via JPL-style EOP files. The state
! (pos & vel) of the center of mass of the Earth with respect to the
! Solar System barycenter is also computed, along with the position of
! the radar station with respect to the Sun. Optionally. the barycentric
! velocities of the station can be computed. The ephemeris time ET
! corresponding to the input UTC time can be either computed or
! specified as an input, according to an input flag.  The expressions
! used to model precession, nutation, and sidereal time can be either
! those of the FK4/B1950 system of constants, or the FK5/ J2000 system,
! as specified in a (mandatory) initialization call to stn_pos_ini.  The
! times utc and et can be either full Julian dates, or offsets from an
! arbitrary reference time t_ref, which is also specified in the call to
! stn_pos_ini.  If the times are full Julian dates, t_ref should be
! specified as zero in the call to stn_pos_ini.
!
!(inertial-frame components of the position and velocity)
! Inputs:
!    utc        - UTC in days relative to t_ref.
!    stn_coords - Geocentric Earth-fixed coordinates of the stations. (AU)
!    radar      - Flag to indicate whether the radar related items will be
!                 computed. Nutation rates are computed only if radar is true.
!    use_eop    - Flag to indicate whether EOP file should be consulted for 
!                 timing and polar motion data. This must be set .true. if 
!                 the radar flag is .true.
!    compute_et - Flag to indicate whether the ephemeris time (ET) is to be 
!                 computed from utc or if it is an input, e.g., for radar 
!                 upleg computation.
!
! Input/Output:
!    et - Ephemeris time in days from t_ref.
!
! Outputs:
!    etmutc     - ET-UTC (sec)
!    x_stn      - Inertial-frame position [and velocity if radar=.true.] of the 
!                 station w.r.t. the Solar-System barycenter at time utc 
!                 (AU, AU/d).
!    x_earth    - Inertial-frame position and velocity of the Earth center 
!                 w.r.t. the Solar-System barycenter at time utc (AU, AU/d).
!    r_stn_sun  - Inertial-frame position of the station w.r.t. the Sun at
!                 time utc (AU).
!    errcod     - Zero indicates success. -1 indicates ETMUT routines were used
!                 instead of EOP.
!
! NOTE: If radar=.false. then the station inertial velocities [x_stn(4:6)] 
!       are not computed.
!$==============================================================================

!^########################## STN_POS_INI #######################################
subroutine stn_pos_ini(j2000, t_ref, errcod)

logical, intent(in) :: j2000
real (kind=wp), intent(in) :: t_ref
integer, intent(out) :: errcod
!===============================================================================
! stn_pos_ini must be called before calling stn_pos. This routine initializes
! the following external routines:
!    eop_mod - loads data and sets t_ref
!    et_minus_ut_mod  - set t_ref
!    etmini - load data file and set internal parameter
!    sdtime - Sidereal Time
!    obliqm - Mean obliquity
!    precmd - Precession
!
! Inputs:
!    j2000 - Flag to indicate FK5/J2000 or FK4/B1950 reference frame
!    t_ref - Reference time, days relative to JD=0.
!
! Outputs:
!    errcod - Zero indicates success.
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../comp/opnetm.f90 <<<<<<<<<<<<<<<<<<<

subroutine opnetm( u )
integer, intent(in) :: u

!-----------------------------------------------------------------------
! Open the ETMUT.dat file and prepare it for reading.
!
! $Id: opnetm.f90,v 1.3 2004/12/02 00:37:15 chesley Exp $
!$----------------------------------------------------------------------
