
>>>>>>>>>>>>>>>>>>> FILE: ../util/astelt.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE ASTELT( ELT, RI, NP, ELA, SIGA)
      DOUBLE PRECISION   ELT(6), RI(*), ELA(4), SIGA(4)
      INTEGER            NP
C
C-----------------------------------------------------------------------
C  ASTELT (ASTeroid ELemenTs) computes the "asteroid" parameters a, P,
C  varpi and M, where "a" is the semi-major axis in AU, "P" is the
C  orbital  period in days,"M" is the mean anomaly, and varpi is the 
C  longitude of perihelion. The routine also computes the sigmas
C  associated with these parameters.
C
C  Inputs:
C   ELT(1..5)          First five of the "comet" elements:
C                       e, eccentricity; 
C                       q, perihelion distance in AU;
C                       Tp, time of perihelion in days wrt to epoch; 
C                       Omega, longitude of ascending node in degrees; 
C                       w, argument of perihelion in degrees.
C   RI(1..NP*(NP+1)/2) Epoch sqrt covariance matrix, upper triangular,
C                      vector stored.
C   NP                 Number of parameters (ie, number of rows and
C                      columns in sqrt covariance matrix).
C
C  Outputs:
C   ELA(1..4)         "Asteroid" elements:
C                       a, semi-major axis in AU; P, orbital period in
C                       days; M, mean anomaly in degrees; varpi, long.
C                       of perihelion in degrees.
C   SIGA(1..4)        Sigmas on the elements in ELA, in the same units 
C                      as ELA.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/bplfrm.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE BPLFRM( V, C)
      DOUBLE PRECISION   V(3), C(3,3)
C
C     SUBROUTINE BPLEPS( CEP, SEP)
      DOUBLE PRECISION   CEP, SEP
C
C-----------------------------------------------------------------------
C  BPLFRM (B-PLane FRaMe) computes the rotation matrix C which rotates
C  the inertial reference frame into the B-plane reference frame.  The
C  B-plane, also known as the target plane, is the plane perpendicular
C  to the incoming asymptote of the hyperbolic path of a spacecraft
C  relative to the target body.  The B-plane frame is defined by the
C  orthonormal set of basis vectors T, R, and S, where
C    T  lies in the B-plane, and is parallel to the ecliptic plane,
C       pointing in the direction of decreasing celestial longitude,
C    R  lies in the B-plane, perpendicular to T and in the direction of
C       decreasing celestial latitude, and
C    S  is directed along the incoming asymptote, perpendicular to the
C       B-plane.
C
C  Input:
C   V(1..3)         Inertial components of a vector along the incoming
C                   asymptote of the spacecraft.  For a massless target
C                   body, the spacecraft velocity vector relative to the
C                   target body at the time of closest approach could be
C                   used.  Any units may be used.
C
C  Output:
C   C( 1..3, 1..3)  Rotation matrix which rotates the inertial frame
C                   into the B-plane frame.
C
C-----------------------------------------------------------------------
C  BPLEPS (B-PLane EPSilon) sets up the sine and cosine of epsilon, the
C  obliquity of the ecliptic relative to the inertial frame.
C
C  Inputs:
C   CEP, SEP        Cosine and sine of obliquity of ecliptic.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/calday.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE CALDAY( YR, MO, DA, TREF, T)
      INTEGER            YR, MO
      DOUBLE PRECISION   DA, TREF, T
C
C-----------------------------------------------------------------------
C  CALDAY (CALendar date to DAY number) converts a specified Gregorian
C  calendar date to a day number relative to a reference Julian day.
C  For example, 1990-Jan-1.0 is day number 14610.0 relative to reference
C  Julian day 2433282.5.  If the reference Julian day is zero, then the
C  day number is the absolute Julian day number.  Relative day numbers
C  allow times to be stored with greater precision.  Negative years are
C  allowed, but they are counted astronomically, so that the year before
C  Year 1 is counted as Year 0, although it is actually the year 1 B.C.
C
C  Inputs:
C   YR    Year in the Gregorian calendar.
C   MO    Month number in the Gregorian calendar.
C   DA    Day number in the Gregorian calendar.  DA is a real variable,
C         and can therefore include a fractional part.
C   TREF  Reference Julian day.
C
C  Output:
C   T     Day number relative to TREF.
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/center.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE CENTER( WIDTH, CSTR)
      INTEGER            WIDTH
      CHARACTER*(*)      CSTR
C
C-----------------------------------------------------------------------
C  CENTER (CENTER) centers the given string within the specified width.
C  If the given string is wider than the specified width, the string is
C  unchanged.  Similarly, if the maximum width of the string is less
C  than the specified width, the string is unchanged.
C
C  Input:
C   WIDTH  The width in which to center the given string.
C
C  Input/Output:
C   CSTR   The string to be centered.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/cetpv.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE CETPV( CE, MU, CEP, SEP, RQPAR, X, ERR)
      DOUBLE PRECISION  CE(6), MU, CEP, SEP, X(3,2,7)
      LOGICAL           RQPAR
      INTEGER           ERR
C
C     SUBROUTINE CETMSG( ERR, MSG)
      CHARACTER*(*)      MSG
C
C-----------------------------------------------------------------------
C  CETPV (Comet Elements to Position and Velocity) converts a specified
C  set of 'comet' orbital elements into components of position and
C  velocity of the orbiting body relative to the center of attraction at
C  current epoch in an inertial frame, given the gravitational parameter
C  MU of the primary body.  The routine also computes the partials of
C  the position and velocity components with respect to the comet
C  elements.  ('Comet' elements use the periapsis distance q and time of
C  periapsis passage Tp instead of semi-major axis a and mean anomaly M).
C
C  Inputs:
C   CE(1..6) Comet orbital elements, defined as follows:
C    CE(1)    eccentricity, e,
C    CE(2)    periapsis distance, q,
C    CE(3)    time of periapsis passage, Tp (relative to current epoch),
C    CE(4)    argument of ascending node, Cap omega (deg),
C    CE(5)    argument of periapsis, small omega (deg),
C    CE(6)    inclination w.r.t. the ecliptic, i (deg).
C   MU       Gravitational parameter of the primary body (MU = GM, where
C            G is gravitational constant and M is mass of the primary).
C   CEP,SEP  Cosine and sine of epsilon, the angle between the ecliptic
C            and the equatorial plane.
C   RQPAR    Set to TRUE to request partials.
C
C   The units of the inputs must be consistent.  For example, if Tp is
C   in days and q is in AU, then MU must be in AU**3/d**2.
C
C  Outputs:
C   X( 1..3, 1, 1)    Position components (same units as q).
C   X( 1..3, 2, 1)    Velocity components (same units as q/Tp).
C   X( 1..3, 1, K+1)  Position partials for Kth parameter.
C   X( 1..3, 2, K+1)  Velocity partials for Kth parameter.
C
C   ERR      Error indicator: nonzero indicates error as follows:
C            1  MU not positive
C            2  e negative
C            3  q negative
C            4  Cannot solve Kepler's equation
C
C  Note:  Kepler's equation is solved for the eccentric anomaly by
C         an iteration which stops when the difference between
C         successive approximations is less than a fixed tolerance
C         or the number of iterations has exceeded a limit.  Both the
C         tolerance and the maximum number of iterations are hard-
C         coded into this routine.  The tolerance is currently set as
C         1E-13 radians and the maximum number of iterations is 100.
C         Different values may be appropriate on other machines.
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/clcase.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE CLCASE( STR)
      CHARACTER*(*)      STR
C
C-----------------------------------------------------------------------
C  CLCASE (Convert to LowerCASE) converts all uppercase letters in a
C  specified character string to lowercase.  This implementation assumes
C  that the character set being used is ASCII.
C
C  Input/Output:
C   STR  Character string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/constants_mod.f90 <<<<<<<<<<<<<<<<<<<

module constants_mod

implicit none
save
private

! Public kind handling variables:
public :: wp,single_precision,double_precision,extended_precision
public :: long_int,short_int

! Public trig constants
public :: pi, two_pi, pi_over_two

! Public unit conversions
public :: deg_per_rad, rad2deg, rad_per_deg, deg2rad, &
          sec_per_rad, rad2sec, rad_per_sec, sec2rad, &
          sec_per_day, day2sec, day_per_sec, sec2day    

! Public astronomical constants
public :: r_earth, au, f_earth, eps_b1950, eps_j2000, &
   gauss_grav_const, gm_sun, v_light, sin_eps, cos_eps, &
   earth_moon_ratio, earth_moon_ratio_p1, au_per_day, gm_conv

!===============================================================================
! constants_mod makes a set of various kind, trig, and astronomical constants 
!   available.
!
!  Modified     Steve Chesley     2001-Nov-07
!     - Added this documentation.
!     - Added pi_over_two.
!  Modified     Steve Chesley     2002-Feb-20
!     - Added earth_moon_ratio, earth_moon_ratio_p1, au_per_day
!  Modified     Steve Chesley     2002-Mar-06
!     - Added gm_conv
!  Modified     Steve Chesley     2002-Mar-28
!     - Added earth_moon_ratio, earth_moon_ratio_p1, au_per_day
!  Modified     Steve Chesley     2003-Sep-17
!     - Tweaked documentation
!  Modified     Steve Chesley     2004-Sep-29
!     - Tweaked documentation
!
! $Id: constants_mod.f90,v 1.7 2004/12/01 23:08:18 chesley Exp $
!===============================================================================

!===============================================================================
!===========                  KIND CONSTANTS                    ================
!===============================================================================

double precision :: test_dp
real :: test_sp

integer, parameter :: single_precision=kind(test_sp)
integer, parameter :: double_precision=kind(test_dp)
integer, parameter :: extended_precision=selected_real_kind(18,200)
integer, parameter :: long_int=selected_int_kind(10)
integer, parameter :: short_int=selected_int_kind(4)
! Assign working precision
integer, parameter :: wp=double_precision
! Test working precision 
! (Will not compile if compiler does not provide desired precision.)
real (kind=wp) :: test_real

!===============================================================================
!===========                  TRIG CONSTANTS                    ================
!===============================================================================

real (kind=wp), parameter :: pi = 3.141592653589793238462643383279502884197_wp
real (kind=wp), parameter :: two_pi = 2.0_wp * pi
real (kind=wp), parameter :: pi_over_two = 0.5_wp * pi

!===============================================================================
!===========              UNIT CONVERSION CONSTANTS             ================
!===============================================================================

! Degrees to radians
real (kind=wp), parameter :: deg_per_rad = 180.0_wp / pi
real (kind=wp), parameter :: rad2deg = deg_per_rad
real (kind=wp), parameter :: rad_per_deg = 1.0_wp / deg_per_rad
real (kind=wp), parameter :: deg2rad  = rad_per_deg

! Arc seconds to radians
real (kind=wp), parameter :: sec_per_rad = 3600.0_wp * deg_per_rad
real (kind=wp), parameter :: rad2sec = sec_per_rad
real (kind=wp), parameter :: rad_per_sec = 1.0_wp / sec_per_rad
real (kind=wp), parameter :: sec2rad  = rad_per_sec

! Time seconds to days
real (kind=wp), parameter :: sec_per_day = 86400.0_wp
real (kind=wp), parameter :: day2sec = sec_per_day
real (kind=wp), parameter :: day_per_sec = 1.0_wp / sec_per_day
real (kind=wp), parameter :: sec2day  = day_per_sec

! big G, GM for misc bodies, au_per_km, rad_earth, v_light, others?

!===============================================================================
!===========                   PLANETARY CONSTANTS              ================
!===============================================================================

! Earth parameters (Geod. ref. sys., 1980) 
real (kind=wp), parameter :: r_earth = 6378.137_wp ! km
real (kind=wp), parameter :: f_earth = 1.0_wp / 298.257222_wp ! flattening

! Astrodynamic Constants not on the Planetary Ephemeris:
! References:
!  1. Explanatory Supplement to the Astronomical Ephemeris,
!       H.M. Nautical Almanac Office, 1961.
!  2. The Astronomical Almanac, 1984, U.S. Government Printing Office.
! Gaussian gravitational constant (Ref 1, p.96)
real (kind=wp), parameter :: gauss_grav_const = .01720209895_wp
! Gravitational parameter (GM) of the Sun, in AU^3/d^2.
real (kind=wp), parameter :: gm_sun = gauss_grav_const ** 2
! Mean obliquities of the ecliptic
! B1950 (Ref 1, p.98)
real (kind=wp), parameter :: time = 0.5_wp ! Private variable
real (kind=wp), parameter :: eps_b1950 = (84428.26_wp - (46.845_wp + &
       (0.0059_wp - 1.81e-3_wp * time) * time) * time) * rad_per_sec
! J2000 (Ref 2, p.S21)
real (kind=wp), parameter :: eps_j2000 = 84381.448_wp * rad_per_sec

!===============================================================================
!===========                   QUASI-CONSTANTS                  ================
!===============================================================================
! These may/must be determined at run time
real (kind=wp) :: au = 1.49597870691e+08_wp ! km
real (kind=wp) :: au_per_day = 1.731456836701389e+03_wp ! (km/sec)/(AU/day)
real (kind=wp) :: gm_conv = 2.229724720949182e-15_wp ! (km^3/sec^2)/(AU^3/day^2)
real (kind=wp) :: earth_moon_ratio = 81.30056_wp
real (kind=wp) :: earth_moon_ratio_p1 = 82.30056_wp
real (kind=wp) :: v_light = 173.1446326846569_wp ! AU/day
real (kind=wp) :: sin_eps, cos_eps

end module constants_mod
!$ End documentation

>>>>>>>>>>>>>>>>>>> FILE: ../util/covtcr.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE COVTCR( COV, SIG, N, COR)
      INTEGER            N
      DOUBLE PRECISION   COV(*), SIG(N), COR(*)
C
C-----------------------------------------------------------------------
C  COVTCR (COVariance To CoRrelation) computes the correlation matrix
C  corresponding to a given covariance matrix.  Both matrices are
C  assumed to be upper triangular vector-stored.  It is also assumed
C  that the standard deviations have already been computed and are
C  available to aid the computation.
C
C  Inputs:
C   COV(1..N*(N+1)/2)  Upper-triangular part of the covariance matrix,
C                      vector-stored.
C   SIG(1..N)          Corresponding sigmas (ie, the square roots of the
C                      diagonal elements of the covariance matrix).
C   N                  Number of rows and columns in COV.
C
C  Output:
C   COR(1..N*(N+1)/2)  Upper-triangular part of the correlation matrix,
C                      vector-stored.  COR may use the same storage area
C                      as input array COV.
C
C  Note:  If any sigma is zero, the corresponding row and column of the
C         correlation matrix is set to zero.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/cucase.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE CUCASE( STR)
      CHARACTER*(*)      STR
C
C-----------------------------------------------------------------------
C  CUCASE (Convert to UpperCASE) converts all lowercase letters in a
C  specified character string to uppercase.  This implementation assumes
C  that the character set being used is ASCII.
C
C  Input/Output:
C   STR  Character string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/eigsym.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE EIGSYM( A, N, NR, W, EVAL, V, ERR)
      INTEGER            N, NR, ERR
      DOUBLE PRECISION   A(*), W(N), EVAL(N), V(NR,N)
C
C-----------------------------------------------------------------------
C  EIGSYM (EIGenvalues and eigenvectors of a real SYMmetric matrix)
C  computes the eigenvalues and eigenvectors of a real symmetric matrix.
C
C  Inputs:
C   A(1..N*(N+1)/2) Upper triangular part of a symmetric matrix of
C                   order N, vector-stored.
C   N               The order of matrix A.
C   NR              Number of rows in V declaration.
C
C  Input/Output:
C   W(1..N)         Work array.
C
C  Outputs:
C   EVAL(1..N)      Eigenvalues of A in ascending order.
C   V(1..NR, 1..N)  Eigenvectors corresponding to the eigenvales
C                   stored in EVAL.
C   ERR             Error indicator: nonzero indicates error as follows:
C                   -1: N <= 0
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/eltj2.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE ELTJ2( ELB5, ELJ2)
C
C     SUBROUTINE ELT50( ELJ2, ELB5)
C
      DOUBLE PRECISION  ELB5(3), ELJ2(3)
C
C-----------------------------------------------------------------------
C  ELTJ2 (ELements To J2000) converts angular orbital elements (longi-
C  tude of ascending node, argument of perihelion, and inclination) from
C  the FK4/B1950 system to the FK5/J2000 system.  The algorithm is due
C  to E.M. Standish.
C
C  Input:
C   ELB5(1..3)  Angular elements wrt mean ecliptic and equinox of B1950
C               on the FK4 system, in the following order:
C       (1)      Longitude of ascending node (deg)
C       (2)      Argument of perihelion (deg)
C       (3)      Inclination (deg)
C
C  Output:
C   ELJ2(1..3)  Angular elements wrt mean ecliptic and equinox of J2000
C               on the FK5 system, in the same order as above.
C
C-----------------------------------------------------------------------
C  ELT50 (ELements To B1950) converts angular orbital elements (longi-
C  tude of ascending node, argument of perihelion, and inclination) from
C  the FK5/J2000 system to the FK4/B1950 system.  The algorithm is due
C  to M.S.W. Keesey.
C
C  Input:
C   ELJ2(1..3)  Angular elements wrt mean ecliptic and equinox of J2000
C               on the FK5 system, in the following order:
C       (1)      Longitude of ascending node (deg)
C       (2)      Argument of perihelion (deg)
C       (3)      Inclination (deg)
C
C  Output:
C   ELB5(1..3)  Angular elements wrt mean ecliptic and equinox of B1950
C               on the FK4 system, in the same order as above.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fcheck.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FCHECK( FNAM, FDESC, RTNNAM)
      CHARACTER*(*)      FNAM, FDESC, RTNNAM
C
C-----------------------------------------------------------------------
C  FCHECK (File CHECK) checks for the existence of a specified file, and
C  if the file does not exist, aborts with an informative message.
C
C  Inputs:
C   FNAM    File name.
C   FDESC   File descriptor, for use in possible abort message.
C   RTNNAM  Name of calling routine, for use in possible abort message.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtbds.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTBDS( STR, LIN, N)
      CHARACTER*(*)      STR, LIN
      INTEGER            N
C
C-----------------------------------------------------------------------
C  FMTSTR (ForMaT Blank Delimited String) formats a specified string as
C  a blank delimited string (eg. putting quotes around it if it contains
C  spaces) and appends it to the character buffer LIN, where on entry, N
C  indicates the number of characters already in LIN, and on exit, N is
C  updated to reflect the appended text.  VARLIST parsing rules are used
C  to determine whether quotes around the string are needed: if the
C  string contains a space, tab, or any of the characters '*', '=', ',',
C  '!', '(', or if its first character is a quote, then the string is
C  enclosed in quotes and quotes within the string are doubled.
C
C  Input:
C   STR   String to be appended to LIN as a blank delimited string.
C
C  Input/Output:
C   LIN   Character buffer to which the blank delimited string is to be
C         appended.
C   N     Number of characters in LIN; updated on output.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtcds.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTCDS( JD, ND, LIN, K)
      DOUBLE PRECISION   JD
      INTEGER            ND, K
      CHARACTER*(*)      LIN
C
C-----------------------------------------------------------------------
C  FMTCDS (ForMaT Calendar Date String) converts a Julian Day number to
C  a calendar date string of the form "yyyy-Mmm-dd.ddd...", where Mmm is
C  a three-character month abbreviation and dd.ddd... is the day of the
C  month, with the specified number of decimal digits, and appends the
C  field to the output line buffer with trailing zeroes on the fraction
C  truncated.  If the number of decimals is specified as -1, the
C  fractional part of the day is simply truncated; otherwise, the field
C  is rounded.
C
C  Inputs:
C   JD     Julian Day number.
C   ND     Number of decimals on day field.
C          If negative, the day field is truncated to integer.
C
C  Input/Outputs:
C   LIN    String containing the result.  LEN(LIN) must be >= K + ND+12
C   K      Number of characters in LIN, updated on output.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtdp.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTDP( VAL, SGNCOL, STR)
      DOUBLE PRECISION  VAL
      LOGICAL           SGNCOL
      CHARACTER*(*)     STR
C
C-----------------------------------------------------------------------
C  FMTDP (ForMaT Double Precision) converts a double precision value to
C  a fixed-length character string using a format which adjusts to
C  provide as many significant digits as possible.  The routine uses
C  either F-format or E-format, whichever allows more significant digits
C  to appear, preferring F-format in case of a tie.  If exponential
C  format is used, one digit appears before the decimal, followed by as
C  many decimal digits as will fit, the letter 'E', and as compact an
C  exponent as possible (eg, E8 or E-6); the decimal point is suppressed
C  if the output string is too short otherwise.  The format adjusts to
C  eliminate leading spaces (except possibly one in the first column, if
C  the calling program reserves column 1 for the sign by setting SGNCOL
C  to TRUE, and the value is non-negative).  The first digit (or the
C  decimal point) appears in column 2 if SGNCOL is TRUE or the value is
C  negative; otherwise (if SGNCOL is FALSE and the value is non-
C  negative), it appears in column 1.  If the value is negative, a minus
C  sign always appears in column 1.  If the output string is too
C  short to accommodate the value under any format, it is filled with
C  asterisks.  If the output string is longer than needed to express the
C  value with maximum precision, the string is padded with trailing
C  blanks.
C
C  Inputs:
C   VAL     Double precision value to be converted.
C   SGNCOL  Set TRUE to reserve the first column for the sign, ie, to
C           force a space in the first column if the value is positive.
C           If SGNCOL is FALSE, digits for non-negative values start in
C           column 1.
C
C  Output:
C   STR     Output string to receive the converted value.
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtdpf.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTDPF( VAL, ND, STR)
      DOUBLE PRECISION   VAL
      INTEGER            ND
      CHARACTER*(*)      STR
C
C-----------------------------------------------------------------------
C  FMTDPF (ForMaT Double Precision, Fixed) converts a double precision
C  value to a fixed-length character string with a specified number of
C  digits right of the decimal point.  If the value overflows, however,
C  the routine uses a format which adjusts to provide as many signifi-
C  cant digits as possible, by using an adjustable F-format or E-format.
C
C  Inputs:
C   VAL   Double precision value to be converted.
C   ND    Number of decimal digits.  If ND = 0, the routine will round
C         to the nearest integer, and use integer format.
C
C  Output:
C   STR     Output string to receive the converted value.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtdpg.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTDPG( VAL, FORM, ND, STR)
      DOUBLE PRECISION   VAL
      INTEGER            ND
      CHARACTER*(*)      FORM, STR
C
C-----------------------------------------------------------------------
C  FMTDPG (ForMaT Double Precision, General) converts a double precision
C  value to one of 2 forms specified with appropriate parameters.
C
C  Inputs:
C   VAL   Double precision value to be converted.
C   FORM  Type of format (with optional sign):
C         "F" Fixed decimal format (F)(+|-)(0)
C         "E" Exponential format (E|e)(+|-)
C   ND    Number of digits after the decimal point for "F", or
C         number of significant digits for "E".
C
C  Output:
C   STR     Output string to receive the converted value.
C           Note: length of STR must be big enough to hold the requested
C                 format.  If not, FMTDPG will abort.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtdps.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTDPS( VAL, NSD, LIN, K)
      DOUBLE PRECISION   VAL
      INTEGER            NSD, K
      CHARACTER*(*)      LIN
C
C-----------------------------------------------------------------------
C  FMTDPS (ForMaT Double Precision, Significant digits specified)
C  converts a double precision value to a compact character string,
C  accurate to the specified number of significant digits, places the
C  result at the current cursor location in the output line, and
C  advances the cursor.  The routine will use exponential format with
C  one digit before the decimal place, or F-format, whichever is more
C  compact.  Trailing zeroes to the right of the decimal point are
C  truncated.
C
C  Inputs:
C   VAL  Double precision value to be formatted.
C   NSD  Number of significant digits.
C
C  Input/Outputs:
C   LIN  String containing the result.
C   K    Number of characters in LIN, updated on output.
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtint.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTINT( NUM, LIN, J)
      INTEGER            NUM, J
      CHARACTER*(*)      LIN
C
C-----------------------------------------------------------------------
C  FMTINT (ForMaT INTeger) converts a specified integer value to a
C  minimal length character string and appends it onto the character
C  string LIN, where on entry, J indicates the number of characters
C  already in LIN.  On exit, J is updated to indicate the new number of
C  characters in LIN.  The number of characters added to the output
C  string varies according to the magnitude of the number, since leading
C  spaces are not stored in the output string; if the value is negative,
C  a leading minus sign is stored in the output string.  This version
C  of FMTINT can handle full 32-bit integers.
C
C  Input:
C   NUM   Integer value to be converted to character.
C
C  Input/Outputs:
C   LIN   Character string in which to place result.
C   J     Number of characters in LIN.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtjds.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTJDS( EPOCH, T, ND, LIN, K)
      DOUBLE PRECISION   EPOCH, T
      INTEGER            ND, K
      CHARACTER*(*)      LIN
C
C-----------------------------------------------------------------------
C  FMTDPS (ForMaT Julian Day, Significance specified) converts a speci-
C  fied time in days past a specified reference Julian Day epoch to a
C  character string accurate to the specified number of decimal digits.
C  The maximum accuracy is preserved by treating the integer and
C  fractional parts of the epoch separately.  Provided the reference
C  epoch is represented exactly internally (ie, it has a simple frac-
C  tional part such as '.5'), this routine can output Julian dates to
C  greater than machine accuracy: eg. 20 digits accuracy on a machine
C  with only 16 digits accuracy.  The field width of the output string
C  is ND + 8.
C
C  Inputs:
C   EPOCH  Epoch time (Julian day).
C   T      Time in days past epoch.
C   ND     Number of decimals.
C
C  Input/Outputs:
C   LIN    String containing the result.
C   K      Number of characters in LIN, updated on output.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtrad.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTRAD( RAD, FMT, ND, RDS)
      DOUBLE PRECISION   RAD
      INTEGER            ND
      CHARACTER*(*)      FMT, RDS
C
C-----------------------------------------------------------------------
C  FMTRAD (ForMaT RADians) converts a specified angle (typically R.A. or
C  declination) into a character string in one of several formats:
C   1.  hours, minutes (or degrees, arcmin)
C   2.  hours, minutes, seconds (or degrees, arcmin, arcsec)
C   3.  degrees
C
C  Inputs:
C   RAD  Angle Value (rad).  Range: [-pi,2*pi] or [-pi,pi]
C   FMT  The format specifier string (R1, R2, R3, D1, D2, or D3
C         where R# implies R.A. and D# implies Dec.).
C   ND   Number of decimal digits.
C             Maximum values of ND for the first 2 formats:
C               FMT   ND
C               R1    6
C               R2    4
C               D1    4
C               D2    3
C
C  Output:
C   RDS  Angle string in one of the following formats:
C           R1. "hh:mm[.m]'"       D1. "+dd mm.[m]'"
C           R2. "hh:mm:ss[.s]"     D2. "+dd mm'ss[.s]""
C           R3. "ddd[.d]"          D3. "+dd[.d]"
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fmtstr.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FMTSTR( STR, LIN, N)
      CHARACTER*(*)      STR, LIN
      INTEGER            N
C
C-----------------------------------------------------------------------
C  FMTSTR (ForMaT STRing) appends a specified string onto the character
C  buffer LIN, where on entry, N indicates the number of characters
C  already in LIN, and on exit, N is incremented by the length of STR.
C
C  Input:
C   STR   Character string to be appended to LIN.
C
C  Input/Output:
C   LIN   Character buffer to which STR is to be appended.
C   N     Number of characters in LIN; updated on output.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/fxtini.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE FXTINI( XSTRT, XSTOP, XSTEP, XTOL, RQMIN, RQMAX,
     .                   ENDPTS, DONE)
      DOUBLE PRECISION   XSTRT, XSTOP, XSTEP, XTOL
      LOGICAL            RQMIN, RQMAX, ENDPTS, DONE
C
C     SUBROUTINE FXTVAL( FVAL, XNXT, DONE)
      DOUBLE PRECISION   FVAL, XNXT
C
C     SUBROUTINE FXTRMA( INDX, FXT, XXT, FND)
      INTEGER            INDX
      DOUBLE PRECISION   FXT, XXT
      LOGICAL            FND
C
C=======================================================================
C  FXTRMA:  This package finds the global extrema (minimum and/or
C  maximum) of a function of one variable f(x) over a specified closed
C  domain [XSTRT,XSTOP].  The user specifies whether or not the
C  endpoints are candidates for these extrema.  Interior extrema are
C  sought using a simple grid of x-values with specified constant
C  stepsize.  The user should specify a stepsize smaller than the
C  expected minimum separation between relative extrema; otherwise the
C  algorithm may not find all the relative extrema.  Each relative
C  extremum is probed to see whether it is a global extremum.  The
C  algorithm used to search for a relative extremum is an adaptation of
C  the golden section search (Ref: Numerical Recipes, 1986, p. 277).
C
C  This package operates in a "handshaking" mode: the user supplies the
C  basic inputs to the initialization routine FXTINI, and then repeated-
C  ly evaluates the function and calls FXTVAL in a DO WHILE (.NOT. DONE)
C  loop.  Thus, the function itself is not passed to this package: the
C  calling program is simply asked to evaluate the function at a series
C  of points.  Each call to FXTVAL supplies one value, and returns the
C  next x-value at which to evaluate the function, or sets DONE= .TRUE.
C  When the algorithm signals that it is done, the results may be
C  obtained by calling FXTRMA.
C
C-----------------------------------------------------------------------
C  FXTINI (Function eXTrema INItialize) initializes this package by
C  specifying the interval over which to probe the function, the step
C  size to use for the probing, the convergence tolerance, the extrema
C  which are sought (minimum and/or maximum), and whether or not the
C  endpoints are candidates for the extrema.  Note that the smallest
C  useful value for XTOL is equal to the square root of machine epsilon
C  times the x-value of the extremum.
C
C  Inputs:
C   XSTRT  Lower bound of x domain and first value at which to evaluate
C          function.
C   XSTOP  Upper bound of x domain.
C   XSTEP  Step size used in searching domain.
C   XTOL   Tolerance in x for relative extrema.
C   RQMIN  Flag set TRUE to request finding minimum.
C   RQMAX  Flag set TRUE to request finding maximum.
C   ENDPTS Flag set TRUE if the endpoints are extrema candidates.
C
C  Output:
C   DONE   Flag set TRUE if the package cannot find extrema because the
C          inputs were incorrect or neither extremum was requested.
C
C-----------------------------------------------------------------------
C  FXTVAL (Function eXTrema VALue) provides a test value for the extrema
C  finding algorithm.  FXTINI must be called first to initialize the
C  package.  The subsequent call to FXTVAL then supplies the value of
C  the function at XSTRT, and returns the next x-value at which to
C  evaluate the function.  FXTVAL should be called repeatedly until
C  the routine signals that it is done.  At this point, FXTRMA can be
C  called to obtain an extremum value.
C
C  Input:
C   FVAL   Function value at current x value.
C
C  Outputs:
C   XNXT   Next x value at which to evaluate the function, provided
C           DONE is false.
C   DONE   Set TRUE if the package has completed finding the extrema.
C
C-----------------------------------------------------------------------
C  FXTRMA (Function eXTReMA) provides the requested extremum value, or
C  indicates that no extremum was found.
C
C  Input:
C   INDX   Index of extremum: 1 for minimum, 2 for maximum.
C
C  Outputs:
C   FXT    Value of requested function extremum.
C   XXT    Abscissa (x-) value for requested function extremum.
C   FND    Set TRUE if the requested extremum was found.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/gettkn.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE GETTKN( STR, K, TYP, TOK, L, KFRST)
      CHARACTER*(*)      STR, TOK
      INTEGER            K, TYP, L, KFRST
C
C-----------------------------------------------------------------------
C  GETTKN (GET TOKeN) scans an input string, starting at a given cursor
C  location, extracts the text for the next unit of data (a "token"),
C  returns an indicator of the type of the token, and leaves the cursor
C  on the next non-blank character following the token.  The calling
C  program can extract all the tokens from a string by calling GETTKN
C  repeatedly.  A token can be any valid FORTRAN numeric field (integer
C  or real), a character string enclosed in quotes, or a name, such as a
C  variable name or body name.  Spaces are not allowed within numeric
C  fields or names. A name token consists of any string of letters and/
C  or digits which is not a numeric field.  A name can also consist of a
C  period followed by an alphanumeric string (eg, .TRUE.): these are
C  called dot names (the optional trailing period is considered to be
C  part of the name).   Except for the periods in dot names, special
C  characters (non-letters and non-digits) are not allowed within names.
C  Lowercase letters in names are converted to uppercase.  Character
C  strings can be delimited by either an apostrophe or double quote.  If
C  it occurs within a string, the delimiting character must be doubled.
C  Note that if the cursor moves off the end of the string at any time,
C  the output cursor value will be larger than the size of the string.
C
C  Inputs:
C   STR   Input character string.
C
C   K     Initial cursor value (ie, character position in STR), K >= 1.
C         GETTKN starts scanning STR at the character pointed to by K.
C         If K > LEN(STR), the routine will return TYP = -1 (see below).
C
C  Outputs:
C   K     Final cursor, set either to the position of the first non-
C         blank character following the token, or, if none was found, a
C         value greater than the length of the input string.  Note that
C         if this character is a delimiter such as a comma, the calling
C         routine should increment K before calling GETTKN again.
C         Otherwise, the next call to GETTKN will return a null name
C         and leave the cursor stuck on the same delimiter.
C
C   TYP   Type of token: TYP = -1: No token found
C                               0: Name
C                               1: Character string
C                               2: Numeric field
C                               3: Integer numeric field
C                               4: Unsigned integer numeric field
C
C   TOK(1:L)  Extracted token text.  If L = 0, TOK is not changed.
C             Lowercase letters in names are converted to uppercase.
C         Dots in dot names are included.  Enclosing delimiters around
C         character strings are stripped; double delimiters within
C         strings are replaced with singles.  Trailing blanks and the
C         delimiter following the token are not copied.
C
C   L     Number of characters in the extracted token.
C
C   KFRST  Cursor pointing to first character of token.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/get_unit_mod.f90 <<<<<<<<<<<<<<<<<<<

module get_unit_mod

implicit none
save
private

! Public procedures
public :: get_unit

integer, parameter :: min_unit=50, max_unit=100

contains

integer function get_unit()

!===============================================================================
! get_unit returns the lowest available unit number in the span 
!   min_unit <= get_unit <= max_unit.
!   The function returns -1 if no unit is that range is unopened.
!-------------------------------------------------------------------------------
! Written      Sometime in 2000         Steve Chesley
!
! $Id: get_unit_mod.f90,v 1.3 2004/12/01 23:08:18 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../util/gram_schmidt_mod.f90 <<<<<<<<<<<<<<<<<<<

module gram_schmidt_mod

use constants_mod

implicit none
save
private

! Public procedures
public :: gram_schmidt

! Public variables
public :: gram_schmidt_msg

! Module variables
character (len=256) :: gram_schmidt_msg

contains

subroutine gram_schmidt(subspace, basis, errcod)

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

!===============================================================================
! gram_schmidt computes an orthonormal basis given vector as one of the 
!     basis directions.
!
! Input:
!   subspace - This is an m x n matrix, the n columns of which span a
!     "constraint" subspace that must be represent by the last n columns 
!     of the new basis.
! Output:
!   basis  - The new basis, given by the columns of the matrix. The last
!      n columns of basis span the subspace given by the input argument.
!   errcod - Zero indicates success.
!
!-------------------------------------------------------------------------------
! Written       Mar 21, 2001         Steve Chesley
! Modified      Mar 12, 2004         Steve Chesley
!    - Major rewrite to use multiple constraint vectors and cope with poor 
!      choice of starting basis vectors.
!
! $Id: gram_schmidt_mod.f90,v 1.3 2004/12/01 22:54:53 chesley Exp $
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../util/hdrlin.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE HDRLIN( PRGVER)
      CHARACTER*(*)      PRGVER
C
C-----------------------------------------------------------------------
C  HDRLIN (HeaDeR LINe) constructs and writes a header line containing
C  the name and version of the program, a mnemonic indicating the
C  operating system and/or compiler, and the run date and time.
C
C  Input:
C   PRGVER  Program name/version string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/hhpost.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE HHPOST( W, IMAXW, IRW, JCW, V, S)
      INTEGER            IMAXW, IRW, JCW
      DOUBLE PRECISION   W(IMAXW,JCW), V(JCW), S(*)
C
C-----------------------------------------------------------------------
C  HHPOST (HouseHolder POST-multiply) triangularizes a rectangular array
C  W by post-multiplying by Householder orthogonal transformations, and
C  places the result in an upper-triangular vector-stored array S.  This
C  version of HHPOST is an improved version of a routine in the ESL
C  library.
C
C  Inputs:
C   IMAXW                 Row dimension of W.
C   IRW                   Number of rows of W to be used in triangulari-
C                         zation (1 <= IRW <= IMAXW).
C   JCW                   Number of columns of W (JCW >= IRW).
C
C  Input/Output:
C   W(1..IMAXW, 1..JCW)   Rectangular sqrt covariance matrix, destroyed
C                         on output.
C
C  Output:
C   V(1..JCW)             Work vector.
C   S(1.. IRW*(IRW+1)/2)  Upper-triangular sqrt covariance matrix,
C                         vector-stored.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/ifnb.f <<<<<<<<<<<<<<<<<<<

      INTEGER FUNCTION IFNB( STR)
      CHARACTER*(*)          STR
C
C-----------------------------------------------------------------------
C  IFNB (Index of First NonBlank) searches a character string and
C  returns the index of the first nonblank character.  If the string
C  is all spaces, the function returns <length of string> + 1.
C
C  Input:
C   STR   Specified string.
C
C  Output:
C   <function value>  Index of first nonblank character, or
C                     LEN(STR) + 1, if STR is all spaces.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/indexx.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE INDEXX(M,N,ARRIN,INDX)

C  This routine returns an index of pointers used to construct a sorted
C  array without swapping variables.  For details, see Press, Flannery,
C  et al on p229+ of 'Numerical Recipes' (1990).
C
C  Inputs:  
C          M       Physical dimension of ARRIN
C          N       Logical dimension of ARRIN
C          ARRIN   The 1-D array to index
C
C  Output:
C          INDX    A 1-D pointer array such that ARRIN(INDX(j)),j=1,n
C                  is in ascending order
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/ischk.f <<<<<<<<<<<<<<<<<<<

      INTEGER FUNCTION ISCHK( ARRAY, N, MNC, KEY)
      INTEGER                 N, MNC
      CHARACTER*(*)           ARRAY(N), KEY
C
C-----------------------------------------------------------------------
C  ISCHC (Index SearCH, Keyword) linearly searches a character string
C  array for the first occurrence of a given keyword value.  The keyword
C  may be abbreviated to a specified minimum of characters.  The search
C  begins with the first element.  The function returns the index of the
C  first matching value, or zero, if the keyword was not found.
C
C  Input:
C   ARRAY(1..N)  Character string array to search.
C   N            Dimension of ARRAY.
C   MNC          Minimum number of characters that must match.
C   KEY          Key value to search for.
C
C  Function returns:
C   Index of first matching value, or 0 if no match occurs.
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/jdcal.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE JDCAL( JD, YR, MO, D)
      DOUBLE PRECISION  JD, D
      INTEGER           YR, MO
C
C-----------------------------------------------------------------------
C  JDCAL (Julian Day to CALendar date) computes the Gregorian calendar
C  date corresponding to a given Julian Day number.  The routine gives
C  correct results for all non-negative Julian Day numbers.  Negative
C  years on output are counted astronomically, so that YR = 0 indicates
C  the year 1 B.C., YR = -1 indicates the year 2 B.C., etc.
C
C  Input:
C   JD     Julian Day number; can include a fractional part of the day.
C
C  Outputs:
C   YR     Year in the Gregorian calendar.
C   MO     Month number in the Gregorian calendar.
C   D      Day number within month, including any fractional part of
C           the day.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/jdcds.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE JDCDS( JD, CDS)
      DOUBLE PRECISION  JD
      CHARACTER*(*)     CDS
C
C-----------------------------------------------------------------------
C  JDCDS (Julian Day to Calendar Date String) converts a Julian Day
C  number to a character string containing the Gregorian calendar date
C  of the form "yyyy Mmm dd.ddd...", where Mmm is a three-character
C  month abbreviation, and dd.ddd... is the day of the month, including
C  fractional part.  The day field contains as many decimal places as
C  will fit in the string.  If the output string is exactly 11
C  characters wide, the fractional part of the day is simply truncated;
C  if it is wider, the field is rounded; if it is shorter, it is filled
C  with asterisks.
C
C  Input:
C   JD   Julian Day number.
C
C  Output:
C   CDS  Calendar date string of the form "yyyy mmm dd.dd...".
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/jdcdts.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE JDCDTS( JD, CDTS)
      DOUBLE PRECISION   JD
      CHARACTER*(*)      CDTS
C
C-----------------------------------------------------------------------
C  JDCDTS (Julian Day to Calendar Date/Time String) converts a Julian
C  Day number to a character string containing the Gregorian calendar
C  date and time in the form "yyyy mmm dd hh:mm:ss.sss...", where yyyy
C  is the year, mmm is a three-character month abbreviation, dd is the
C  day of the month, hh is the hour, mm is the minute, and ss.sss... is
C  the second, including the fractional part.  The seconds field
C  contains as many decimal places as will fit in the string.  If the
C  output string is exactly 20 characters long, the time is rounded to
C  the nearest second; if the output string is 17-19 characters long,
C  the time is rounded to the nearest minute; if it is shorter than 17
C  characters, it is filled with asterisks.
C
C  Input:
C   JD    Julian Day number.
C
C  Output:
C   CDTS  Calendar date/time string of the form
C          "yyyy mmm dd hh:mm:ss.sss...".
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/jdtfmt.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE JDTFMT( JD, DFMT, CHRMO, SEP, TFMT, STR)
      DOUBLE PRECISION   JD
      INTEGER            DFMT, TFMT
      LOGICAL            CHRMO
      CHARACTER          SEP*1, STR*(*)
C
C-----------------------------------------------------------------------
C  JDTFMT (Julian Day date/Time ForMaT) formats a Julian day number into
C  a date/time character string in a user-specified format.  The output
C  string has the general format <DATE><SEP><TIME>, where any of these
C  fields may be null.  <SEP> is the separation character, which is
C  included only if both <DATE> and <TIME> are non-null.  The format of
C  the date is determined by the format code DFMT, which indicates the
C  desired order of subfields year, month, and day, and flag CHRMO,
C  which indicates whether the month should be formated as 3 characters
C  (MMM) or 2 digits (mm), as follows:
C
C               CHRMO = .TRUE.      CHRMO = .FALSE.
C               --------------      ---------------
C   DFMT = 0     <null>               <null>
C          1     yyyy-MMM-dd          yyyy-mm-dd
C          2     MMM dd, yyyy         mm/dd/yyyy
C          3     dd-MMM-yyyy          dd-mm-yyyy
C
C  If the time field is to be omitted (TFMT = 0), and the day subfield
C  is rightmost (DFMT = 1), the fraction of the day will be included
C  in the date field, with digits continuing to the end of the specified
C  string.
C
C  The format of the time is determined by format code TFMT, which
C  indicates whether the time should be just hours and fractions
C  thereof, hours and minutes and fractions of minutes, or hours,
C  minutes and seconds and fractions of seconds, as follows:
C
C   TFMT = 0     <null>
C          1     hh.hhh...
C          2     hh:mm.mmm...
C          3     hh:mm:ss.sss...
C
C  Digits in the time field continue to the end of the specified string.
C
C  Inputs:
C   JD        Julian day, including fraction of a day, if desired.
C   DFMT      Date format:  0: no date field;  1: year, month, day;
C               2: month, day, year;  3: day, month, year.
C   CHRMO     Flag: TRUE to request character month,
C                   FALSE for numeric month.
C   SEP       Separation character.
C   TFMT      Time format:  0: no time field;  1: hours only;
C               2: hour and minutes only;  3: hours, minutes, seconds.
C
C  Output:
C   STR       Formatted date/time string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/julian.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE JULIAN( YR, MO, D, JD)
      INTEGER            YR, MO
      DOUBLE PRECISION   D, JD
C
C-----------------------------------------------------------------------
C  JULIAN computes the Julian Day number corresponding to a given
C  calendar date in the Gregorian calendar.  Negative years are allowed,
C  but they are counted astronomically, so that the year before Year 1
C  is counted as Year 0, although it is actually the year 1 B.C.
C
C  Inputs:
C   YR     Year in the Gregorian calendar.
C   MO     Month number in the Gregorian calendar.
C   D      Day number within month, including a possible fractional
C           part.
C
C  Output:
C   JD     Julian Day number.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/linini.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE LININI( UNIT, LINSIZ)
      INTEGER            UNIT, LINSIZ
C
C     SUBROUTINE LINCHK( LIN, LLEN)
      CHARACTER*(*)      LIN
      INTEGER            LLEN
C
C     SUBROUTINE LINWRT( LIN, LLEN)
C
C     SUBROUTINE LINMRG( MARGIN)
      INTEGER            MARGIN
C
C=======================================================================
C  This is a package of routines for managing the packing of variable-
C  length fields into a line-oriented output file.  The calling routine
C  is expected to maintain the line buffer and current line length (ie,
C  it must append the new field onto the line buffer and increment the
C  current line length).  After it adds a new field to the line buffer,
C  the calling routine should call LINCHK to see whether the line has
C  overflowed the desired line length.  If this has occurred, LINCHK
C  outputs the first part of the line buffer, up to but not including
C  the new field, and then shifts the new field to the left margin of
C  the line buffer.
C
C-----------------------------------------------------------------------
C  LININI (LINe INItialize) initializes this package with the unit
C  number for outputting lines and the desired line size.  The left
C  margin is initialized to zero.  Note that the desired line size is
C  smaller than the size of the line buffer.  Note also that some output
C  lines might exceed the desired line size, since each output line
C  contains at least one field, and the sum of the margin size plus
C  field size may exceed the desired line size.  LIMINI must be called
C  before any other routines in this package.  The output file need not
C  be open when LININI is called, but should be open before LIMADD is
C  called.
C
C  Inputs:
C   UNIT    Unit number for output file.
C   LINSIZ  Desired line size, in characters, including the left margin.
C
C-----------------------------------------------------------------------
C  LINCHK (LINe CHecK) checks the line buffer after a field has been
C  added to it.  If the new length of the buffer exceeds the desired
C  line length for the file, the contents of the line buffer up to but
C  not including the new field are output to the file, and the line
C  buffer is reset to remove the fields which were written.
C
C  Input/Output:
C   LIN     Line buffer.
C   LLEN    Current line length.
C
C-----------------------------------------------------------------------
C  LINWRT (LINe WRiTe) terminates the current line after the field added
C  in the last call to LINCHK, writes the line to the output file, and
C  re-initializes the line and line length according to the current
C  margin setting.
C
C  Input/Output:
C   LIN     Line buffer, current on input, re-initialized on output.
C
C  Output:
C   LLEN    Line length of re-initialized line.
C
C-----------------------------------------------------------------------
C  LINMRG (LINe MaRGin) sets the margin to be used the next time the
C  line buffer is reset or re-initialized (ie, it does not affect the
C  current contents of the buffer).  The default margin is zero.
C
C  Input:
C   MARGIN  New margin setting.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/lubksb.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE LUBKSB( A, N, NP, INDX, B)
      INTEGER            N, NP, INDX(N)
      DOUBLE PRECISION   A(NP,N), B(N)
C
C-----------------------------------------------------------------------
C  LUBKSB (LU BacK SuBstitution) performs the back substitution step in
C  the solution of the linear equation Ax=b, where A has already been
C  expanded into its LU decomposition via routine LUDCMP, and b is a
C  specified column matrix.
C
C  Inputs:
C   A(1..N, 1..N)  LU decomposition of the original A matrix.
C   N              Size of A.
C   NP             Row dimension of A.
C   INDX           Record of row permutations.
C
C  Input/Output:
C   B              Right-hand side of Ax=b on input,
C                  solution vector x on output.
C
C  Ref.: Numerical Recipes, p.36
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/ludcmp.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE LUDCMP( A, N, NP, INDX, SGN)
      INTEGER            N, NP, INDX(N)
      DOUBLE PRECISION   A(NP,N), SGN
C
C-----------------------------------------------------------------------
C  LUDCMP (LU DeCoMPosition) performs an LU decomposition on a specified
C  N x N matrix A, and replaces A with the result.  INDX is a record of
C  the row permutations effected by the partial pivoting.  SGN is the
C  parity of the row permutations unless A is singular, in which case
C  SGN is zero.  This routine is used in combination with LUBKSB to
C  solve linear equations or to invert a matrix.
C
C  Input/Output:
C   A(1..N, 1..N)  Input square matrix to be decomposed on input;
C                  LU decomposition of A on output.
C
C  Inputs:
C   N              Size of A.
C   NP             Row dimension of A.
C
C  Outputs:
C   INDX(1..N)     Record of row permutations.
C   SGN            Parity of row permutations (+1. or -1.), used in
C                   computation of the determinant, or
C                  zero if A is singular.
C
C  Ref.: Numerical Recipes, p.35
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/magnit.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE MAGNIT( X, XH, RHO, R, MP, COMET, VM, FND)
      DOUBLE PRECISION   X(6), XH(6), RHO, R, MP(5), VM
      LOGICAL            COMET, FND
C
C     SUBROUTINE MAGINI( VMAGIN)
      DOUBLE PRECISION   VMAGIN
C
C-----------------------------------------------------------------------
C  MAGNIT (MAGNITude) computes the V-band magnitude given heliocentric
C  and topocentric states and distances, magnitude parameters, and a
C  flag set true if the object is a comet.  MAGINI must be called first.
C
C  Inputs:
C   X(1..6)   Topocentric state.
C   XH(1..6)  Heliocentric state.
C   RHO       Topocentric distance (AU).
C   R         Heliocentric distance (AU).
C   MP(1..5)  Magnitude parameters:
C               1: M1 (or H if not a comet)
C               2: M2 (or G if not a comet)
C               3: K1
C               4: K2
C               5: Phase coefficient
C   COMET     Flag set true if the object is a comet (i.e. uses comet
C             nuclear magnitude law).
C
C  Outputs:
C   VM        V-magnitude estimate (1 x 10^9 if not available).
C   FND       True if not fainter than specified limit (see MAGINI).
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/mapsrc.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE MAPSRC( RI, NP, UV, NC, W)
      INTEGER            NP, NC
      DOUBLE PRECISION   RI(*), UV(6,NC), W(6,NP)
C
C-----------------------------------------------------------------------
C  MAPSRC (MAP Square Root Covariance) maps a square root covariance
C  matrix from one epoch and/or set of orbital elements to another,
C  given the mapping matrix UV.  It is assumed that only the first 6
C  elements are affected by this mapping.  The output mapped sqrt
C  covariance matrix is left rectangular (ie, is not triangularized).
C
C  Inputs:
C   RI(1..NP*(NP+1)/2) Input sqrt covariance matrix, upper triangular,
C                      vector-stored.
C   NP                 Number of parameters (ie, number of rows and
C                      columns in sqrt covariance matrix).
C   UV( 1..6, 1..NC)   Mapping matrix.
C   NC                 Number of columns in mapping matrix; NC <= NP.
C
C  Output:
C   W( 1..6, 1..NP)    First 6 rows of mapped sqrt covariance matrix.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/matvec_util_mod.f90 <<<<<<<<<<<<<<<<<<<

module matvec_util_mod

use constants_mod

implicit none
save
private ! by default

! Public procedures:
public :: mat_inv, vec_mag, vec_unit, vec_cross, vsmat2mat, vsmat_vec_mul, &
   vsmat_mat_mul, vsmat_diagmat_mul

!===============================================================================
! A mishmash of convenience routines for dealing with vectors and matrices. 
!   Not all of these are very efficient so use with caution.
!
! mat_inv           - computes the inverse of a square matrix
! vec_mag           - computes the magnitude of a vector
! vec_cross         - computes the cross product of two three-vectors
! vec_unit          - computes a unit vector
! vsmat2mat         - converts a vector-stored matrix to a matrix
! vsmat_vec_mul     - vector-stored matrix & vector multiplication
! vsmat_mat_mul     - vector-stored matrix & matrix multiplication
! vsmat_diagmat_mul - vector-stored matrix & diagonal matrix multiplication
!
! NOTE: Error handling is a bit irregular. Instead of setting a public msg 
!   variable when the error flag is set, a warning message is output to 
!   STDOUT and STDERR.
!
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
! Modified      many times from 2000-2004
! Modified      Nov 30, 2004       Steve Chesley
!     - Improved documentation
!-------------------------------------------------------------------------------
! $Id: matvec_util_mod.f90,v 1.9 2004/12/01 23:07:40 chesley Exp $
!$==============================================================================

!^########################## MAT_INV ###########################################
subroutine mat_inv(mat, inv, error)

real (kind=wp), intent(in)  :: mat(:,:)
real (kind=wp), intent(out) :: inv(:,:)
integer,        intent(out) :: error

!===============================================================================
! mat_inv computes the inverse of a square matrix using LU decomposition.
!-------------------------------------------------------------------------------
! Written       Aug 17, 2001       Steve Chesley
!$==============================================================================

!^########################## VEC_MAG ###########################################
function vec_mag(vec)

real (kind=wp), intent(in) :: vec(:)
real (kind=wp)             :: vec_mag

!===============================================================================
! vec_mag computes the magnitude of a vector
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
!$==============================================================================

!^########################## VEC_CROSS #########################################
function vec_cross(vec1, vec2)

real (kind=wp), intent(in) :: vec1(:), vec2(:)
real (kind=wp)             :: vec_cross(3)

!===============================================================================
! vec_mag computes the cross product of two three-vectors. Array entries 
!   above three are ignored. If the input arrays are shorter than three 
!   then the results will be unreliable.
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
!$==============================================================================

!^########################## VEC_UNIT ##########################################
function vec_unit(vec)

real (kind=wp), intent(in) :: vec(:)
real (kind=wp)             :: vec_unit(size(vec))

!===============================================================================
! vec_unit computes a unit vector
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
!$==============================================================================

!^########################## VSMAT2MAT #########################################
function vsmat2mat(vsmat, n, symmetric)

real (kind=wp), intent(in)    :: vsmat(:)
integer,        intent(in)    :: n
logical, optional, intent(in) :: symmetric
real (kind=wp)                :: vsmat2mat(n,n)

!===============================================================================
! vsmat2mat (vector-stored matrix to matrix) returns the vsmat as an n x n
!    matrix. If symmetric is .true. then the lower triangle is initialized 
!    symmetrically, otherwise the lower triangle elements are zero. If the 
!    vsmat vector is longer than necessary, the extra elements at the end are 
!    ignored.
!
! Inputs: 
!    vsmat - Vector-stored matrix. The vsmat elements are assumed to be stored 
!        in the following scheme:
!             1  2  4  7  11 ...
!                3  5  8  12 ...
!                   6  9  13 ...
!                     10  14 ...
!                         15 ...
!                            ...
!    n - the dimension of vsmat
! Optional Input:
!    symmetric - Indicates whether vsmat is symmetric or upper triangular.
!       (Default = .false.)
!-------------------------------------------------------------------------------
! Written       Dec 14, 2000       Steve Chesley
!$==============================================================================

!^########################## VSMAT_VEC_MUL #####################################
subroutine vsmat_vec_mul(vsmat, vec, order, vec_out, error, utm)

real (kind=wp),    intent(in)  :: vsmat(:)
real (kind=wp),    intent(in)  :: vec(:)
character (len=*), intent(in)  :: order
real (kind=wp),    intent(out) :: vec_out(:)
integer,           intent(out) :: error
logical, optional, intent(in)  :: utm

!===============================================================================
! vsmat_vec_mul (vector-stored matrix & vector multiplication) 
!
! Inputs: 
!    vsmat -   Vector-stored matrix.
!    vec -     Vector to be multiplied
!    order -   'PRE' to premultiply (vec*vs), 'POST' to postmultiply (vs*vec).
! Output:
!    vec_out - The product of the multiplicands, we hope.
!    error -   Zero indicates succes.
! Optional Input:
!    utm -     .true. to indicate that vsmat is an Upper Triangular Matrix, 
!              .false. if vsmat is a symmetric matrix. (Default is .true.)
!-------------------------------------------------------------------------------
! Written       Aug 29, 2001       Steve Chesley
!$==============================================================================

!^########################## VSMAT_MAT_MUL #####################################
subroutine vsmat_mat_mul(vsmat, mat, order, mat_out, error)

real (kind=wp),    intent(in)  :: vsmat(:)
real (kind=wp),    intent(in)  :: mat(:,:)
character (len=*), intent(in)  :: order
real (kind=wp),    intent(out) :: mat_out(:,:)
integer,           intent(out) :: error

!===============================================================================
! vsmat_mat_mul (vector-stored matrix & matrix multiplication) 
!
! Inputs: 
!    vsmat - Vector-stored matrix.
!    mat - Matrix to be multiplied
!    order - 'PRE' to premultiply (mat*vs), 'POST' to postmultiply (vs*mat).
! Output:
!    mat_out - The product of the multiplicands, we hope.
!    error - zero indicates succes.
!-------------------------------------------------------------------------------
! Written       Aug 29, 2001       Steve Chesley
!$==============================================================================

!^########################## VSMAT_DIAGMAT_MUL #################################
subroutine vsmat_diagmat_mul(vsmat, diagmat, order, vsmat_out, error)

real (kind=wp),    intent(in)  :: vsmat(:)
real (kind=wp),    intent(in)  :: diagmat(:)
character (len=*), intent(in)  :: order
real (kind=wp),    intent(out) :: vsmat_out(:)
integer,           intent(out) :: error

!===============================================================================
! vsmat_diagmat_mul (vector-stored matrix & diagonal matrix multiplication) 
!
! Inputs: 
!    vsmat - Vector-stored matrix.
!    diagmat - Vector representing diagonal matrix to be multiplied
!    order - 'PRE' to premultiply (mat*vs), 'POST' to postmultiply (vs*mat).
! Output:
!    vsmat_out - The product of the multiplicands, we hope.
!    error - zero indicates succes.
!-------------------------------------------------------------------------------
! Written       Aug 29, 2001       Steve Chesley
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../util/moid_mod.f90 <<<<<<<<<<<<<<<<<<<

module moid_mod

use constants_mod
use matvec_util_mod

implicit none
save
private ! by default

! Public procedures:
public :: comp_moid, get_extrema

!===============================================================================
! The moid_mod uses the Sitarski method (Acta Astronomica v. 18, p. 171 1968) 
! to compute the extrema of the relative distances of two orbits. 
!
! The idea is to call comp_moid with two sets of elements, returning the MOID, 
! or the number of extrema obtained. Given the number of extrema the details 
! for each extremum can be obtained with a subsequent call to get_extrema.
!
! $Id: moid_mod.f90,v 1.6 2004/12/01 22:57:25 chesley Exp $
!
! Modified    2004 Dec 01   Steve Chesley
!   - Cleanup and documentation. Removed get_dist routine.
!
!$==============================================================================

!^########################## COMP_MOID #########################################
subroutine comp_moid(elements1, elements2, moid, true_anom1, true_anom2, &
   n_extrema, n_step)

real (kind=wp),           intent(in)  :: elements1(6)
real (kind=wp),           intent(in)  :: elements2(6)
real (kind=wp), optional, intent(out) :: moid
real (kind=wp), optional, intent(out) :: true_anom1
real (kind=wp), optional, intent(out) :: true_anom2
integer,        optional, intent(out) :: n_extrema
integer,        optional, intent(in)  :: n_step

!===============================================================================
! comp_moid computes the extrema in distance between two orbits, returning 
!   any of several optional outputs
!
! Inputs: 
!   elements1  - First set of elements in the following order
!                  [e q (AU), tp (not used), Om (deg), w (deg), i (deg)]
!   elements2  - Second set of elements, same as elements1
!
! Optional Outputs:
!   moid       - MOID (AU)
!   true_anom1 - True anomaly of MOID point for 1st orbit
!   true_anom2 - True anomaly of MOID point for 2nd orbit
!   n_extrema  - Number of extrema found. These can be retrieved with a 
!                subsequent call to get_extrema.
! Optional Inputs:
!   n_step     - The number of steps to take in the one-dimensional search.
!                Default is 1000.
!$==============================================================================

!^########################## GET_EXTREMA #######################################
subroutine get_extrema(type_vec, dist_vec, v1_vec, v2_vec)

character (len=3),        intent(out)  :: type_vec(:)
real (kind=wp),           intent(out)  :: dist_vec(:)
real (kind=wp),           intent(out)  :: v1_vec(:)
real (kind=wp),           intent(out)  :: v2_vec(:)

!===============================================================================
! get_extrema returns arrays containing the tabulated extrema information 
!   saved on the previous call to comp_moid.
!
! Outputs:
!   type_vec - Type of extremum: MIN, MAX, or SAD
!   dist_vec - Distances of the extrema
!   v1_vec   - Orbit 1 true anomalies
!   v2_vec   - Orbit 2 true anomalies
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../util/obscode_mod.f90 <<<<<<<<<<<<<<<<<<<

module obscode_mod

use constants_mod
use get_unit_mod

implicit none
save
private ! by default

! Public procedures:
public :: obscode_ini, obscode_index, observatory_ini

! Public data types:
public :: observatory_type

! Public parameters:
public :: len_obs_name, len_coords_str, max_code

! Public variables:
public :: obscode_msg, sites, num_obscode

! Field sizes
integer, parameter :: len_obs_name       = 60
integer, parameter :: len_coords_str     = 15

! The maximum number of observatories that can be handled
integer, parameter :: max_code = 2000

! A defined data type:
type observatory_type
  character (len=3)              :: code
  character (len=3)              :: fmt_type
  character (len=len_coords_str) :: coords_str(3)
  real (kind=wp)                 :: coords_xyz(3)
  character (len=len_obs_name)   :: name
end type observatory_type

! Retrievable error message
character (len=180) :: obscode_msg

! Array of observatory data, and the number of loaded observatories
type (observatory_type) :: sites(max_code)
integer :: num_obscode

!===============================================================================
! The obscode_mod module and associated routines do the following tasks:
!      - the module defines the observatory_type data type
!      - obscode_ini loads the observatory information into the 'sites' 
!        observatory_type array 
!      - obscode_index returns the 'sites' index location of a given obscode
!      - observatory_ini intializes the internal numerical values for a 
!        given observatory_type variable
!
!-------------------------------------------------------------------------------
! Written       May 19, 2000       Steve Chesley
! Modified      Jul 31, 2001       Steve Chesley
!      - fixed coords_str format statements
! Modified      Jul 31, 2001       Steve Chesley
!      - documentation
!
! $Id: obscode_mod.f90,v 1.6 2004/12/01 23:06:46 chesley Exp $
!$==============================================================================

!^############# OBSCODE_INI ####################################################
subroutine obscode_ini(errcod)

integer, intent(out) :: errcod

!===============================================================================
! obscode_ini loads the observatory coordinates into memory. The routine
!   needs two files: 'obscode.mpc', which is the text portion of the MPC
!   Obscode.html, and 'obscode.jpl', which is the JPL file (radar). If a
!   code is in both files the JPL file will overwrite the MPC data. The
!   routine looks for the files in $CAET_DATA if defined, else in $OBSDIR
!   if defined, else in the working directory.
!
! Output:
!   errcod   - error code. Zero indicates success.
!$==============================================================================

!^########### OBSCODE_INDEX ###################################################
integer function obscode_index(obscode)

character (len=*), intent(in) :: obscode

!===============================================================================
! obscode_index returns the index of a particular array in the sites 
!   array. Returns zero if obscode is not found.
!
! Input:
!   obscode - the three character observatory code.
!$==============================================================================

!^########### OBSERVATORY_INI ##################################################
subroutine observatory_ini(observatory, errcod)

type (observatory_type), intent(inout) :: observatory
integer,                 intent(out)   :: errcod

!===============================================================================
! observatory_ini fills in the coords_xyz field in an observatory_type variable
!
! Input/Output:
!   observatory - observatory data type variable that needs to be initialized
! Output:
!   errcod      - error code. Zero indicates success.
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../util/obstj2.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE OBSTJ2( TO, RA50, DC50, RAJ2, DCJ2)
C
C     SUBROUTINE OBST50( TO, RAJ2, DCJ2, RA50, DC50)
      DOUBLE PRECISION   TO, RAJ2, DCJ2, RA50, DC50
C
C
C-----------------------------------------------------------------------
C  OBSTJ2 (OBServation To J2000) converts an astrometric observation
C  from the FK4/B1950 system to the FK5/J2000 system.  Elliptic aberra-
C  tions (E-terms) are assumed to have been removed from the input.
C  The output R.A. and Dec. may overwrite the input R.A. and Dec.
C
C  Inputs:
C   TO    Time of the observation (JED).
C   RA50  Right ascension in the FK4/B1950 system (rad).
C   DC50  Declination in the FK4/B1950 system (rad).
C
C  Outputs:
C   RAJ2  Right ascension in the FK5/J2000 system (rad); range [-pi,pi].
C   DCJ2  Declination in the FK5/J2000 system (rad).
C
C-----------------------------------------------------------------------
C  OBST50 (OBServation To B1950) converts an astrometric observation
C  from the FK5/J2000 system to the FK4/B1950 system.  Elliptic aberra-
C  tions (E-terms) are not included in the output FK4/B1950 position.
C  The output R.A. and Dec. may overwrite the input R.A. and Dec.
C
C  Inputs:
C   TO    Time of the observation (JED).
C   RAJ2  Right ascension in the FK5/J2000 system (rad).
C   DCJ2  Declination in the FK5/J2000 system (rad).
C
C  Outputs:
C   RA50  Right ascension in the FK4/B1950 system (rad); range [-pi,pi].
C   DC50  Declination in the FK4/B1950 system (rad).
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/orbunc.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE ORBUNC( P, EC, SIGP, SIGTP, RUNOFF, U)
      DOUBLE PRECISION   P, EC, SIGP, SIGTP, RUNOFF
      INTEGER            U
C
C-----------------------------------------------------------------------
C  ORBUNC (ORBit UNCertainty) computes the Minor Planet Center
C  uncertainty parameter U given an object's period and eccentricity, as
C  well as the uncertainty in period and perihelion time.
C
C  Inputs:
C   P       Orbital period (d).
C   EC      Orbital eccentricity.
C   SIGP    One-sigma uncertainty in P (d).
C   SIGTP   One-sigma uncertainty in time of perihelion (d).
C
C  Output:
C   RUNOFF  In-orbit longitude runoff (arc-sec/decade).
C   U       Uncertainty parameter.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/parsob.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PARSOB( OBJ, NAME, IAUNUM, DESIG)
      CHARACTER*(*)      OBJ, NAME, DESIG
      INTEGER            IAUNUM
C
C-----------------------------------------------------------------------
C  PARSOB (PARSe OBject) parses a full object name field, as given in an
C  orb file, and returns its components: the IAU number, name, and, if
C  present, the  designation.
C
C  Input:
C   OBJ     Object name field, as given in the orb file.
C
C  Output:
C   NAME    IAU name, without the IAU number or designation (except for
C           unnumbered comets, for which a designation in parentheses
C           is not stripped from the name because it may be required to
C           establish uniqueness).
C           Example: 'Ida', 'P/Halley', '1986 JK'
C   IAUNUM  IAU number of object.
C           If the object is unnumbered, a value of zero is returned.
C   DESIG   Designation, if present; otherwise blanks.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/posfrm.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE POSFRM( R, RE, C)
      DOUBLE PRECISION   R(3), RE(3), C(3,3)
C
C-----------------------------------------------------------------------
C  POSFRM (Plane-Of-Sky FRaMe) computes the rotation matrix C which
C  rotates the inertial reference frame into the plane-of-sky reference
C  frame as viewed from the Earth.  This frame is defined by the
C  orthonormal set of basis vectors A, D, and R, where
C    A  lies in the plane of sky as viewed from the Earth, pointing in
C       the direction of increasing right ascension,
C    D  also lies in the plane of sky as viewed from the Earth, perpen-
C       dicular to A, in the direction of increasing declination, and
C    R  is directed along the range vector from the Earth to the
C       orbiting body, perpendicular to the plane of sky.
C
C  Input:
C   R(1..3)         Inertial-frame position of the orbiting body w.r.t.
C                   the Sun.  Any units may be used.
C   RE(1..3)        Inertial-frame position of the Earth w.r.t. the
C                   Sun.  Must use the same units as R.
C
C  Output:
C   C( 1..3, 1..3)  Rotation matrix which rotates the inertial frame
C                   into the plane-of-sky frame.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/prtels.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PRTELS( SYS, SOLN, EPOCH, X, DX, PNM, NP, NAMDX, NELA,
     .                   ELA, DXA)
C
      USE CONSTANTS_MOD, ONLY: GM_SUN, GM_CONV
C
      INTEGER            NP, NELA
      CHARACTER*(*)      SYS, SOLN, NAMDX, PNM(NP)*3
      DOUBLE PRECISION   EPOCH, X(NP), DX(NP), ELA(4), DXA(4)
C
C-----------------------------------------------------------------------
C  PRTELS (PRinT ELementS) outputs a columnar summary of the orbital
C  elements and other estimated parameters in an easily readable format,
C  along with a second column indicating either the accuracy of the
C  elements or the amount by which they have changed from a previous
C  solution.  Each element is output in its own appropriate format.  For
C  example, the time of periapsis passage is printed in both Julian day
C  and calendar date formats.  The elements are output in the order in
C  which they are stored, which can be any order.  The routine
C  distinguishes each element by its name, which is passed through
C  argument PNM.  Semi-major axis A and mean anomaly MA will also be
C  printed, if they are passed to this routine via AMA(1) and AMA(2),
C  respectively.  They are not printed if AMA(1)=0.
C
C  Inputs:
C   SYS         Reference system name for elements; eg. 'J2000'.
C   SOLN        Solution name.
C   EPOCH       Epoch of the elements (Julian day, ET).
C   X(1..NP)    Parameter values, including the orbital elements.
C               Angular elements (eg. inclination) are assumed to be in
C               degrees.
C   DX(1..NP)   Values for the second column (eg., corrections or
C               standard deviations).  Values for angular parameters are
C               assumed to be in radians and are converted to degrees
C               when output.
C   PNM(1..NP)  Names of the parameters.
C   NP          Number of parameters.
C   NAMDX       Heading for the second column.
C   NELA        No. of elements in ELA; supported values: 0, 3, or 4.
C   ELA(1..NELA) Auxiliary elements, semi-major axis a, period P, 
C                Long. perihelion w~, and mean anomaly M.  List may be 
C                truncated.
C   DXA(1..NELA) Auxiliary element sigmas, in same units as ELA.  List 
C                may be truncated.
C
C  Output:
C   Summary on unit 6.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/prtmap.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PRTMAP( TMAP, W, RI, NP, PNM, PNX, MFRAME, SAVMAP)
      INTEGER            NP
      DOUBLE PRECISION   TMAP, W(6,NP), RI(*)
      CHARACTER*3        PNM(NP), PNX(6), MFRAME*(*)
      LOGICAL            SAVMAP
C
C-----------------------------------------------------------------------
C  PRTMAP (PRinT MAP output) computes and outputs the covariance matrix
C  associated with a given square root covariance matrix which has been
C  mapped to the map time and rotated into a particular frame.  The
C  routine also outputs the corresponding correlation matrix, the
C  standard deviations, and the ellipse parameters (semi-major axis,
C  semi-minor axis, and orientation angle theta).  The name of the
C  mapping frame is output in the heading above each of these
C  quantities, and the rows and columns are labelled with the
C  appropriate component names.  If requested, this routine also writes
C  the mapped covariance matrix and associated information to a file
C  in NAMELIST format.
C
C  Inputs:
C   TMAP                Map time (Julian day, ET).
C   W( 1..6, 1..NP)     First 6 rows of the mapped and rotated sqrt
C                       covariance matrix (km & km/s).
C   RI(1..NP*(NP+1)/2)  The original sqrt covariance matrix at epoch
C                       from which the sqrt covariance of non-state
C                       parameters is obtained.  This is upper-
C                       triangular vector-stored.
C   NP             Number of parameters.
C   PNM(1..NP)     Parameter names of epoch state parameters.
C   PNX(1..6)      Parameter names for mapping reference frame.
C   MFRAME         Name of mapping reference frame.
C   SAVMAP         Set TRUE to request that the mapped covariance and
C                  associated parameters be written to a file in
C                  NAMELIST format.
C
C  Outputs:
C   1. Output on unit number 6.
C   2. If SAVMAP is set, a file containing the mapped covariance and
C      associated parameters.
C
C  Note:  The error ellipse is formed by projecting the error ellipsoid
C         defined by the top left 3x3 partition of the covariance matrix
C         onto the plane spanned by the first two axes of the mapping
C         frame.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/prtrow.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PRTROW( U, N, NAM, TTL)
      INTEGER            N
      DOUBLE PRECISION   U(N)
      CHARACTER*3        NAM(N), TTL*(*)
C
C-----------------------------------------------------------------------
C  PRTROW (PRinT ROW) outputs a vector horizontally across a line,
C  complete with column labels and a title.  Up to nine columns are
C  output on each line.  If more are needed, the vector is output in
C  segments of up to nine columns each.
C
C  Inputs:
C   U(1..N)    Vector of values to be output.
C   N          Number of components in U.
C   NAM(1..N)  Parameter names of the components of U.
C   TTL        Title to be output above vector.
C
C  Output:
C   Output on unit 6.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/prtsta.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PRTSTA( U, N, NAM, TTL)
      INTEGER            N
      DOUBLE PRECISION   U(N)
      CHARACTER*3        NAM(N), TTL*(*)
C
C-----------------------------------------------------------------------
C  PRTSTA (PRinT STA) outputs a vector horizontally across a line,
C  complete with column labels and a title.  Up to nine columns are
C  output on each line.  If more are needed, the vector is output in
C  segments of up to nine columns each.
C
C  Inputs:
C   U(1..N)    Vector of values to be output.
C   N          Number of components in U.
C   NAM(1..N)  Parameter names of the components of U.
C   TTL        Title to be output above vector.
C
C  Output:
C   Output on unit 6.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/prtutm.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PRTUTM( U, N, NAM, TTL)
      DOUBLE PRECISION   U(*)
      INTEGER            N
      CHARACTER*3        NAM(N), TTL*(*)
C
C-----------------------------------------------------------------------
C  PRTUTM (PRinT Upper Triangular Matrix) outputs an upper triangular,
C  vector-stored matrix, complete with column and row labels and title.
C  Up to nine columns are output across each line.  If the matrix
C  requires more than nine columns, it is output in vertical segments of
C  up to nine columns each.
C
C  Inputs:
C   U(1..N*(N+1)/2)   An upper triangular vector-stored matrix.
C   N                 Number of rows and columns in matrix U.
C   NAM(1..N)         Parameter names associated with the rows and
C                     columns of U.
C   TTL               Title to be output above matrix.
C
C  Output:
C   Output on unit 6.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/pvtce.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PVTCE( X, MU, CEP, SEP, CE, ERR)
      DOUBLE PRECISION  X(6), MU, CEP, SEP, CE(6)
      INTEGER           ERR
C
C     SUBROUTINE PVTMSG( ERR, MSG)
      CHARACTER*(*)      MSG
C
C-----------------------------------------------------------------------
C  PVTCE (Position and Velocity To Comet Elements) converts the position
C  and velocity of an orbiting body at the current epoch, in equatorial
C  inertial-frame components, to the ecliptic 'comet' orbital elements
C  for that body, given the gravitational parameter MU of the primary
C  body. ('Comet' elements are those which use the periapsis distance q
C  and time of periapsis Tp instead of the semi-major axis a and mean
C  anomaly M.)
C
C  Inputs:
C   X(1..3)  Equatorial inertial-frame position components at current
C            epoch.
C   X(4..6)  Equatorial inertial-frame velocity components at current
C            epoch.
C   MU       Gravitational parameter of the primary body (MU = GM, where
C            G is gravitational constant and M is mass of the primary).
C   CEP,SEP  Cosine and sine of epsilon, the angle between the
C            ecliptic and the equatorial plane.
C
C   The units of the inputs must be consistent.  For example, if
C   X(1..3) is in AU and X(4..6) in AU/d, MU must be in AU**3/d**2.
C
C  Outputs:
C   CE(1..6) Comet orbital elements, defined as follows:
C    CE(1)    eccentricity, e,
C    CE(2)    periapsis distance, q (same units as X(1..3)),
C    CE(3)    time of periapsis passage, Tp (relative to current epoch,
C              same units as X(1..3)/X(4..6)),
C    CE(4)    argument of ascending node, Cap omega (deg),
C    CE(5)    argument of periapsis, small omega (deg), and
C    CE(6)    inclination w.r.t. the ecliptic, i (deg).
C
C   ERR      Error indicator: nonzero indicates error as follows:
C            1  MU is not positive
C            2  Rectilinear (straight-line) orbit
C
C  Notes:  1. If the orbit lies in the ecliptic plane (i=0 or 180 deg),
C             the argument of the ascending node is set to zero.
C          2. If the orbit is circular (e=0), the current position is
C             used as periapsis.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/pvtsbe.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE PVTSBE( X, GM, RP, SBE, ERR)
      DOUBLE PRECISION   X(6), GM, RP, SBE(6)
      INTEGER            ERR
C
C     SUBROUTINE PVTSBP( PSB)
      DOUBLE PRECISION   PSB(6,6)
C
C     SUBROUTINE PVTSBS( SVC)
      DOUBLE PRECISION   SVC(3)
C
C-----------------------------------------------------------------------
C  PVTSBE (Position and Velocity To Scaled B-plane Elements) converts
C  the state of an orbiting body (components of position and velocity)
C  relative to a primary body into scaled B-plane elements relative to
C  the primary.  The orbiting body must be on a hyperbolic trajectory
C  relative to the primary.  Any units may be used for the inputs, as
C  long as they are consistent; output units will then be consistent
C  with the input units.
C
C  Inputs:
C   X(1..6)  State of orbiting body relative to primary: three position
C             components in length units followed by three velocity
C             components in velocity units.
C
C   GM       Gravitational parameter of primary body (length units times
C             velocity units squared).  GM = G*M, where G is the
C             gravitational constant and M is the mass of the primary).
C   RP       Scaling radius (length units).
C
C
C  Outputs:
C   SBE(1..6) Scaled B-plane elements, defined as follows:
C    SBE(1)    b_pr, r-component of scaled b-vector (length units)
C    SBE(2)    b_pt, t-component of scaled b-vector (length units)
C    SBE(3)    Linearized time of flight, T_L (time units relative to
C               current epoch, where time units are length units over
C               velocity units).
C    SBE(4)    r-component of s-vector
C    SBE(5)    t-component of s-vector
C    SBE(6)    C3 = square of v-infinity (velocity units squared)
C
C   ERR      Error indicator: nonzero value indicates error, as follows:
C            1  MU not positive
C            2  Body at origin (r = 0)
C            3  Orbit not hyperbolic (1/a >= 0)
C            4  Asymptote along z-axis (b-plane not defined)
C
C-----------------------------------------------------------------------
C  PVTSBP (PVTSBe Partials) computes the partials of the scaled b-plane
C  elements with respect to the state for the elements computed in the
C  previous call to PVTSBE.
C
C  Output:
C   PSB(1..6,1..6)  Partials of scaled b-plane elements (rows) with
C                   respect to state (columns).
C
C-----------------------------------------------------------------------
C  PVTSBS (PVTSBe S-vector) returns the components of the s-vector
C  computed in the previous call to PVTSBE.  The reference frame is the
C  same as used for X in the call to PVTSBE.
C
C  Output:
C   SVC(1..3)  S-vector components.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/razel0.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE RAZEL0( TIST, LAST, LASTD, PHI, AZ, EL)
      DOUBLE PRECISION   TIST(6), LAST, LASTD, PHI, AZ(2), EL(2)
C
C---------------------------------------------------------------------
C RAZEL0 - Compute azimuth, elevations and their rates. This routine
C generally follows Escobal's notation and procedure; p 405 thru 407 
C of "Methods of Orbit Determination".
C
C Inputs:
C   TIST  -- Topocentric (apparent) state vector
C   LAST  -- Local Apparent Sidereal Time, radians
C   LASTD -- Sidereal rate, rad/d 
C   PHI   -- Geodetic latitude of observer, radians
C  
C Outputs:
C   AZ    -- AZ(1) = Azimuth, North->East, radians
C            AZ(2) = Azimuthal rate, rad/d
C   EL    -- EL(1) = Elevation above local horizon, radians
C            EL(2) = Elevation rate, rad/d
C
C Author: Jon Giorgini (1996)
C Minor modifications by Alan Chamberlin (September 1996)
C$--------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/razel.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE RAZEL( TIST, LAST, LASTD, PHI, AZ, EL)
      DOUBLE PRECISION  TIST(6), LAST, LASTD, PHI, AZ(2), EL(3)
C
C-----------------------------------------------------------------------
C RAZEL - Compute azimuth, elevations and their rates.
C
C Input:
C   TIST   Topocentric (apparent) state vector
C   LAST   Local Apparent Sidereal Time, radians
C   LASTD  Sidereal rate, rad/s 
C   PHI    Geodetic latitude of observer, radians
C  
C Output:
C   AZ  AZ(1)= Azimuth, North->East, radians
C       AZ(2)= Azimuthal rate, radians/sec
C   EL  EL(1)= Elevation above local horizon, radians
C       EL(2)= Elevation rate, rad/sec
C       EL(3)= Refraction corrected elevation, radians
C
C Refraction references:
C    [1] J. Meeus, "Astronomical Algorithms", 1991, p. 101-102
C    [2] T. Saemundsson, Sky and Telescope, July, 1986, p.70
C
C Author: J. Giorgini (1996)                                  Horizons
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/rdapr.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE RDAPR( U, R, NP, PNM)
      INTEGER           U, NP
      DOUBLE PRECISION  R(*)
      CHARACTER         PNM(NP)*3
C
C-----------------------------------------------------------------------
C  RDAPR (ReaD A PRiori) reads an a-priori square root information array
C  which has been saved in a file in type 66 format by the routine
C  WRSRI.  The sqrt information array is comprised of the sqrt
C  information matrix and the right-hand side of the a-priori data
C  equation.  The number of parameters and the parameter names are also
C  read from the file.  The file must already be open for input, but is
C  closed by this routine. The file is written in ASCII not binary.
C
C  Input:
C   U                  Fortran unit number to use for input.
C
C  Outputs:
C   R(1..NP*(NP+3)/2)  A-priori sqrt information array.
C   NP                 Number of parameters in a-priori sqrt information
C                      array.
C   PNM(1..NP)         Parameter names for a-priori information.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/reord.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE REORD( RO, NAMO, LO, R, NAM, L, A, IR)
      DOUBLE PRECISION  RO(*), R(*), A(IR,*)
      INTEGER           LO, L, IR
      CHARACTER*(*)     NAMO(LO), NAM(L)
C
C-----------------------------------------------------------------------
C  REORD (REORDer) reorders the columns of sqrt information matrix RO
C  with parameter names NAMO according to a new list of parameters NAM.
C  The result is then repacked into matrix R.
C
C  Inputs:
C   RO(*)        Sqrt info matrix, old order.
C   NAMO(1..LO)  Parameter names, old order.
C   LO           Number of parameters old parameter list.
C   NAM(1..L)    New list of parameters.
C   L            Number of parameters in new list.
C   IR           Row dimension of work area A.
C
C  Outputs:
C   R(*)           Sqrt info matrix, new order.
C   A(1..IR,1..L)  Work area.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/rotsrc.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE ROTSRC( WI, NP, C, W)
      INTEGER            NP
      DOUBLE PRECISION   WI(6,NP), C(3,3), W(6,NP)
C
C-----------------------------------------------------------------------
C  ROTSRC (ROTate Square Root Covariance) rotates a square root
C  covariance matrix from one frame to another, given the rotation
C  matrix C which rotates the first frame to the second.  The first 6
C  parameters of the covaraince are assumed to be the position and
C  velocity components in the input frame.  The square root covariance
C  matrix is assumed to be partitioned with the first 6 rows stored
C  separately in a rectangular array WI.  Only these first 6 rows are
C  affected by the change of frames.
C
C  Inputs:
c   WI( 1..6, 1..NP)  First 6 rows of the input square root covariance
C                     matrix.
C   NP                Number of parameters (ie, number of columns in the
C                     sqrt covariance matrix).
C   C( 1..3, 1..3)    Rotation matrix which rotates the input frame into
C                     the output frame.
C
C  Output:
C   W( 1..6, 1..NP)   First 6 rows of the output sqrt covariance matrix.
C                     W must NOT occupy the same storage area as WI.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/rottj2.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE ROTTJ2( V50, VJ2)
      DOUBLE PRECISION   V50(3), VJ2(3)
C
C     SUBROUTINE ROTT50( VJ2, V50)
C
C-----------------------------------------------------------------------
C  ROTTJ2 (ROTate To J2000) rotates a vector from the FK4/B1950 system
C  to the FK5/J2000 system.  The output vector may overwrite the input.
C
C  Inputs:
C   V50(1..3)  FK4/B1950 components of vector.
C
C  Ouptut:
C   VJ2(1..3)  FK5/J2000 components of vector; may overwrite V50(*).
C
C-----------------------------------------------------------------------
C  ROTT50 (ROTate To B1950) rotates a vector from the FK5/J2000 system
C  to the FK4/B1950 system.  The output vector may overwrite the input.
C
C  Inputs:
C   VJ2(1..3)  FK5/J2000 components of vector.
C
C  Ouptut:
C   V50(1..3)  FK4/B1950 components of vector; may overwrite VJ2(*).
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/rtnfrm.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE RTNFRM( X, C)
      DOUBLE PRECISION   X(6), C(3,3)
C
C-----------------------------------------------------------------------
C  RTNFRM (RTN FRaMe) computes the rotation matrix C which rotates the
C  inertial reference frame into the RTN reference frame.  The RTN frame
C  is defined by the orthonormal set of basis vectors R, T, and N, where
C   R  is directed along the radius vector from the Sun to the orbiting
C      body,
C   T  lies in the orbital plane, perpendicular to R and in the general
C      direction of the velocity (ie, the transverse direction), and
C   N  is directed normal to the orbital plane with a sense established
C      by the right-hand rule (with fingers in the direction of motion).
C
C  Input:
C   X(1..6)         Inertial-frame components of position and velocity:
C                      X(1..3) are the components of position, and
C                      X(4..6) are the components of velocity.
C                   Any units may be used.
C
C  Output:
C   C( 1..3, 1..3)  Rotation matrix which rotates the inertial frame
C                   into the RTN frame.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/s3par.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE S3PAR( X, MU, T, PX, ERR)
      DOUBLE PRECISION  X(6), MU, T, PX(6,6)
      INTEGER           ERR
C
C-----------------------------------------------------------------------
C  S3PAR (Set III PARameters) computes the partial derivatives of
C  position and velocity of an orbiting body at the current time with
C  respect to Set III variables at epoch.
C
C  Inputs:
C   X(1..6)   State (position and velocity) at the current time,
C             expressed in the inertial frame.
C   MU        Gravitational parameter of the primary body.
C   T         Time since epoch.
C
C   The units of the inputs must be consistent: if T is in days, and
C   X(1..3) is in AU, then X(4..6) must be in AU/d and MU must be in
C   AU^3/d^2.
C
C  Outputs:
C   PX(1..6, 1..6)  Partial derivatives of X with respect to Set III
C                   variables (DMW, DP, DQ, EDW, DA, DE) at epoch.
C   ERR       Error indicator: a nonzero value indicates an error:
C              ERR = 1  Rectilinear (straight-line) orbit,
C              ERR = 2  Circular orbit,
C              ERR = 3  Parabolic orbit.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/srccov.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE SRCCOV( W, RI, NP, P, SIG)
      INTEGER            NP
      DOUBLE PRECISION   W(6,NP), RI(*), P(*), SIG(NP)
C
C-----------------------------------------------------------------------
C  SRCCOV (Square Root Covariance to COVariance) computes the covariance
C  matrix corresponding to a given square root covariance matrix which
C  has been partitioned so that the first 6 rows are stored in a
C  rectangular array W, and the remaining rows are obtained from the
C  corresponding rows of an array RI, which is upper-triangular and
C  vector-stored.  The routine places the upper triangular part of the
C  covariance into a vector-stored array P, and the sigmas (square roots
C  of the diagonal elements) into SIG.
C
C  Inputs:
C   W( 1..6, 1..NP)     First 6 rows of the sqrt covariance matrix.
C   RI(1..NP*(NP+1)/2)  A full upper-triangular vector-stored matrix
C                       from which rows 7..NP of the sqrt covariance
C                       matrix are obtained.
C   NP                  Number of parameters (ie, number of rows and
C                       columns in the sqrt covariance matrix).
C
C  Outputs:
C   P(1..NP*(NP+1)/2)   Upper triangular part of the covariance matrix,
C                       vector-stored.
C   SIG(1..NP)          Corresponding standard deviations.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/tbstat.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE TBSTAT( CE, EPOCH, X, XT, RHO)
      DOUBLE PRECISION   CE(6), EPOCH, X(6), XT(6), RHO
C
C     SUBROUTINE TBSINI( TOBS, TREF, XOBS, J2000, CIN)
      DOUBLE PRECISION   TOBS, TREF, XOBS(6), CIN
      LOGICAL            J2000
C
C-----------------------------------------------------------------------
C  TBSINI (Two-Body State INItialize).
C
C  Inputs:
C   TOBS        Time of observation (ET) relative to TREF.
C   TREF        Reference JD.
C   XOBS(1..6)  Heliocentric-equatorial observer state.
C   J2000       System to use for equinox and ecliptic:
C                .TRUE.  = J2000 
C                .FALSE. = B1950
C   CIN         Speed of light (AU/d).
C
C-----------------------------------------------------------------------
C  TBSTAT (Two-Body STATe).
C
C  Input:
C   CE(1..6)  Comet-style orbital elements.
C   EPOCH     JD of the element's epoch.
C
C  Output:
C   X(1..6)   Heliocentric-equatorial state of object referred to J2000
C             (or B1950) and approximately corrected for light-time.
C   XT(1..6)  Topocentric-equatorial state of object (X - XO).
C   RHO       Topocentric distance (AU).
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/timfmt.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE TIMFMT( DT, S, DFMT, CHRMO, SEP, TFMT, STR)
      DOUBLE PRECISION   S
      INTEGER            DT(5), DFMT, TFMT
      CHARACTER          SEP*1, STR*(*)
      LOGICAL            CHRMO
C
C-----------------------------------------------------------------------
C  TIMFMT (TIMe ForMaT) formats a "standard" Data/Time array into a
C  string in a format determined by user inputs.  The output string has
C  the general format <DATE><SEP><TIME>, where any of these fields may
C  be null.  <SEP> is the separation character, which is included only
C  if both <DATE> and <TIME> are non-null.  The format of the date is
C  determined by the format code DFMT, which indicates the desired
C  order of subfields year, month, and day, and flag CHRMO, which
C  indicates whether the month should be formated as 3 characters (MMM)
C  or 2 digits (mm), as follows:
C
C               CHRMO = .TRUE.      CHRMO = .FALSE.
C               --------------      ---------------
C   DFMT = 0     <null>               <null>
C          1     yyyy-MMM-dd          yyyy-mm-dd
C          2     MMM dd, yyyy         mm/dd/yyyy
C          3     dd-MMM-yyyy          dd-mm-yyyy
C          4     yyyy MMM dd          yyyy mm dd
C
C  If the time field is to be omitted (TFMT = 0), and the day subfield
C  is rightmost (DFMT = 1), the fraction of the day will be included
C  in the date field, with digits continuing to the end of the specified
C  string.
C
C  The format of the time is determined by format code TFMT, which
C  indicates whether the time should be just hours and fractions
C  thereof, hours and minutes and fractions of minutes, or hours,
C  minutes and seconds and fractions of seconds, as follows:
C
C   TFMT = 0     <null>
C          1     hh.hhh...
C          2     hh:mm.mmm...
C          3     hh:mm:ss.sss...
C
C  Digits in the time field continue to the end of the specified string.
C
C  Inputs:
C   DT(1..5)  Date/Time array in "standard" order: yr, mo, day, hr, min.
C   S         Seconds.
C   DFMT      Date format:  0: no date field;  1: year, month, day;
C               2: month, day, year;  3: day, month, year;  4: same as
C               (1) but the dashes are spaces.
C   CHRMO     Flag: TRUE to request character month,
C                   FALSE for numeric month.
C   SEP       Separation character.
C   TFMT      Time format:  0: no time field;  1: hours only;
C               2: hour and minutes only;  3: hours, minutes, seconds.
C
C  Output:
C   STR       Formatted date/time string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/trhhp.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE TRHHP( MXNR, FROW, LROW, W, RI)
      INTEGER           MXNR, FROW, LROW
      DOUBLE PRECISION  W(MXNR,*), RI(*)
C
C-----------------------------------------------------------------------
C  TRHHP (TRiangularize via HouseHolder transformations Post-multiplied)
C  triangularizes a square root covariance data array consisting of a
C  square matrix W and an upper triangular portion RI.  The upper
C  triangular result overwrites RI.  Householder orthonormal transform-
C  ations are used to perform the triangularization, and special effort
C  is made to avoid overflow.  The argument FROW indicates the first row
C  and column of RI to be triangularized; columns of RI before column
C  FROW are not affected by the triangularization.  Row 1 of W aligns
C  with row FROW of RI.  The argument LROW indicates the last row of RI
C  to be triangularized; rows and columns of RI beyond LROW are not
C  affected by the triangularization.  On input, the upper triangular
C  diagonal partition of RI from rows FROW through LROW is assumed to
C  be zero; this portion of RI is filled in during the triangularization,
C  and rows 1..FROW-1 above it are modified appropriately.
C
C  Inputs:
C   MXNR                    Row dimension of W.
C
C   FROW                    First row of RI to be triangularized;
C                           (1 <= FROW <= LROW).
C
C   LROW                    Last row of RI to be triangularized.
C                           (LROW-FROW+1 <= MXNR).
C
C  Input/Output:
C   W(1..NB,1..NB)          Square data array, where NB = LROW-FROW+1.
C                           Row 1 of W corresponds to row FROW of RI.
C
C   RI(1..LROW*(LROW+1)/2)  Upper triangular square root covariance
C                           matrix, vector stored.  On input, the upper
C                           triangular diagonal partition from rows FROW
C                           to LROW is assumed to be zero, but the rows
C                           above it (1..FROW-1) are used.  On output,
C                           only columns FROW..LROW of RI are changed:
C                           the rest of RI is unaffected by the triangu-
C                           larization.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/wldairm.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE WLDAIRM( RELV, AIRMASS )
C
C-----------------------------------------------------------------------
C WLDAIRM -- Returns approximate relative optical airmass given apparent
C refracted elevation angle. Uses the approximation function of [1]. 
C
C Inputs:
C   RELV    -- Apparent refracted elevation in degrees
C
C Output:
C   AIRMASS -- Approximate relative airmass
C
C Reference:
C
C   [1] Kasten, F., Young, A., "Revised Optical Air Mass Tables and
C       Approximation Formula", Applied Optics, vol 28, no. 22, 
C       p. 4735-4738, Nov. 15, 1989.
C
C Author: Jon Giorgini (1996)
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/wrcov.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE WRCOV( TMAP, P, NP, PNAM)
      DOUBLE PRECISION  TMAP, P(*)
      INTEGER           NP
      CHARACTER*3       PNAM(NP)
C
C-----------------------------------------------------------------------
C  WRCOV (WRite COVariance) writes a file containing the mapped
C  covariance matrix P as well as the map time, the number of
C  parameters, and the names of the parameters, all in NAMELIST format.
C  This data can then be used by other programs.
C
C  Inputs:
C   TMAP               Map time (Julian Day, ET).
C   P(1..NP*(NP+1)/2)  Upper-triangular part of the mapped covariance
C                      matrix, vector-stored.
C   NP                 Number of parameters.
C   PNAM(1..NP)        Parameter names.
C
C  Output:
C   File containing the mapped covariance matrix and associated
C   information.
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/wrsri.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE WRSRI( U, R, NP, PNM)
      INTEGER           U, NP
      DOUBLE PRECISION  R(*)
      CHARACTER*3       PNM(NP)
C
C-----------------------------------------------------------------------
C  WRSRI (WRite Square Root Information) writes the square root
C  information array to a file so that it may be used as the a-priori in
C  another run.  The sqrt information array includes the right-hand side
C  of the data equation.  The number of parameters and parameter names
C  are also saved in the file.  The file must already be open for
C  output, but is closed by this routine.
C
C  Inputs:
C   U                  Fortran unit number to use for output.
C   R(1..NP*(NP+3)/2)  Square root information array.
C   NP                 Number of parameters.
C   PNM(1..NP)         Parameter names.
C
C  Output:
C   A file containing the number of parameters, the parameter names, and
C   the square root information array.  The file is written in type 66
C   format, but in ASCII not binary.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/xtrdbl.f <<<<<<<<<<<<<<<<<<<

      DOUBLE PRECISION FUNCTION XTRDBL( STR)
      CHARACTER*(*)                     STR
C
C-----------------------------------------------------------------------
C  XTRDBL (eXTRact DouBLe) extracts a double precision value from a
C  character string and returns the value as the value of the function.
C  Leading spaces and an optional sign are permitted.  However, the
C  value must be contained within the first 99 characters of the string.
C  The string is assumed to be a valid FORTRAN double precision numeric
C  string.  Note that this function may not be referenced in a WRITE
C  statement since it is implemented using an internal READ statement.
C
C  Inputs:
C   STR   Input character string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/xtrint.f <<<<<<<<<<<<<<<<<<<

      INTEGER FUNCTION XTRINT( STR)
      CHARACTER*(*)            STR
C
C-----------------------------------------------------------------------
C  XTRINT (eXTRact INTeger) extracts an integer value from a character
C  string and returns the value as the value of the function.  Leading
C  spaces and an optional sign are permitted.  Extraction continues
C  until either a non-digit delimiter is found or the end of the string
C  is reached.  This routine is specialized for machines with 32-bit
C  two's complement integers.  The maximum absolute value that can be
C  extracted is 2^31-1 (2147483647).  If no value is found, XTRINT
C  returns a value of -2^31 (-2147483648).
C
C  Inputs:
C   STR   Input character string.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../util/xztopo.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE XZTINI( RE, F)
      DOUBLE PRECISION   RE, F
C
C     SUBROUTINE XZTOPO( X, Z, LATG, HT)
      DOUBLE PRECISION   X, Z, LATG, HT
C
C-----------------------------------------------------------------------
C  XZTINI (X-Z to Topo INItialize) initializes the XZTOPO package.
C
C  Inputs:
C   RE    Equatorial radius of ellipsoid (km).
C   F     Flattening of ellipsoid:  f = 1 - sqrt( 1 - e^2 ),
C          where e is the eccentricity of the ellipsoid cross-section.
C
C-----------------------------------------------------------------------
C  XZTOPO (X-Z coordinates to TOPO coordinates) converts a Cartesian
C  position above an oblate planet into geodetic latitude and height
C  above the reference ellipsoid.
C
C  Reference:
C   Escobal, P.R., Methods of Orbit Determination, 1965, p.398.
C
C  Inputs:
C   X     Cartesian coordinate parallel to equator (km); X >= 0.
C   Z     Cartesian coordinate parallel to polar axis (km).
C
C  Outputs:
C   LATG  Geodetic latitude (rad).
C   HT    Height above reference ellipsoid (km).
C
C$----------------------------------------------------------------------
