sun267.tex

\documentclass[twoside,11pt,nolof]{starlink}

% ? Specify used packages
% ? End of specify used packages

% -----------------------------------------------------------------------------
% ? Document identification
% Fixed part
\stardoccategory  {Starlink User Note}
\stardocinitials  {SUN}
\stardocsource    {sun\stardocnumber}
\stardoccopyright{%
Copyright \copyright\ 2012 Science and Technology Facilities
  Council.\\ Copyright \copyright\ 2014 Cornell University.\\
Copyright \copyright\ 2015 Tim Jenness}

% Variable part - replace [xxx] as appropriate.
\stardocnumber    {267.3}
\stardocauthors   {Tim Jenness}
\stardocdate      {2025 April 17}
\stardoctitle     {PAL --- Positional Astronomy Library}
\stardocversion   {0.9.10}
\stardocmanual    {Programmer's Manual}
\stardocabstract  {%
PAL provides a subset of the Fortran SLALIB library but written in C
using the SLALIB C API. Where possible the PAL routines are
implemented using the C SOFA/ERFA library. It is provided with a GPLv3 license.
}
% ? End of document identification
% -----------------------------------------------------------------------------

\begin{document}
\scfrontmatter

% ? Main text

\section{Introduction}

This library provides a C library designed as a API-compatible
replacement for the C SLALIB library (SUN/67) and uses a GPL licence so is
freely redistributable. Where possible the functions call equivalent
SOFA routines (Hohenkerk, C., 2011, Scholarpedia, \textbf{6}, \emph{11404})\footnote{or equivalent ERFA routines.}
and use current IAU 2006 standards. This means that any
functions that rely on nutation or precession will return slightly
different answers to the SLA functions.

\section{Citing PAL}

If you use PAL in your work please consider citing it. The description paper
for PAL is: \emph{PAL: A Positional Astronomy Library}, Jenness, T. \& Berry, D. S.,
in \emph{Astronomical Data Anaysis Software and Systems XXII}, Friedel, D. N. (ed),
ASP Conf.\ Ser. \textbf{475}, p307.

\clearpage
\appendix
\section{\label{APP:SPEC}Function Descriptions}

By default PAL is set up to use the ERFA variant of SOFA. ERFA is an
approved redistribution of the SOFA code using a BSD-license and
renamed function calls. Whereas SOFA routines have a \texttt{iau} prefix
the ERFA equivalents have a \texttt{era} prefix. The PAL build script
will try to detect which of ERFA and SOFA is available. Wherever SOFA
is mentioned in this document the ERFA equivalent can be substituted.

ERFA can be obtained from \htmladdnormallink{https://github.com/liberfa/erfa}.

\subsection{SOFA Mappings}

The following table lists PAL/SLA functions that have direct
replacements in SOFA. Whilst these routines are implemented in the PAL
library using SOFA new code should probably call SOFA directly.

\begin{tabbing}
\hspace*{2cm}\=\hspace*{3cm}\= \kill
SLA/PAL \>  SOFA \\
\texttt{palCldj} \> \texttt{iauCal2jd} \\
\texttt{palDbear} \> \texttt{iauPas} \\
\texttt{palDaf2r} \> \texttt{iauAf2a} \\
\texttt{palDav2m} \>  \texttt{iauRv2m} \\
\texttt{palDcc2s} \>  \texttt{iauC2s} \\
\texttt{palDcs2c} \> \texttt{iauS2c} \\
\texttt{palDd2tf} \> \texttt{iauD2tf}\\
\texttt{palDimxv} \> \texttt{iauTrxp}\\
\texttt{palDm2av} \> \texttt{iauRm2v}\\
\texttt{palDjcl} \> \texttt{iauJd2cal}\\
\texttt{palDmxm} \> \texttt{iauRxr}\\
\texttt{palDmxv} \> \texttt{iauRxp}\\
\texttt{palDpav} \> \texttt{iauPap}\\
\texttt{palDr2af} \> \texttt{iauA2af}\\
\texttt{palDr2tf} \> \texttt{iauA2tf}\\
\texttt{palDranrm} \> \texttt{iauAnp}\\
\texttt{palDsep} \> \texttt{iauSeps}\\
\texttt{palDsepv} \> \texttt{iauSepp}\\
\texttt{palDtf2d} \> \texttt{iauTf2d}\\
\texttt{palDtf2r} \> \texttt{iauTf2a}\\
\texttt{palDvdv} \> \texttt{iauPdp}\\
\texttt{palDvn} \> \texttt{iauPn}\\
\texttt{palDvxv} \> \texttt{iauPxp}\\
\texttt{palEpb} \> \texttt{iauEpb}\\
\texttt{palEpb2d} \> \texttt{iauEpb2d}\\
\texttt{palEpj} \> \texttt{iauEpj}\\
\texttt{palEpj2d} \> \texttt{iauEpj2jd}\\
\texttt{palEqeqx} \> \texttt{iauEe06a}\\
\texttt{palFk5hz} \> \texttt{iauFk5hz} \textit{also calls iauEpj2jd}\\
\texttt{palGmst} \> \texttt{iauGmst06}\\
\texttt{palGmsta} \> \texttt{iauGmst06}\\
\texttt{palHfk5z} \> \texttt{iauHfk5z} \textit{also calls iauEpj2jd}\\
\texttt{palRefcoq} \> \texttt{iauRefco}\\
\end{tabbing}

%% include palsun.tex

\sstroutine{
   palAddet
}{
   Add the E-terms to a pre IAU 1976 mean place
}{
   \sstdescription{
      Add the E-terms (elliptic component of annual aberration)
      to a pre IAU 1976 mean place to conform to the old
      catalogue convention.
   }
   \sstinvocation{
      void palAddet ( double rm, double dm, double eq,
                      double $*$rc, double $*$dc );
   }
   \sstarguments{
      \sstsubsection{
         rm = double (Given)
      }{
         RA without E-terms (radians)
      }
      \sstsubsection{
         dm = double (Given)
      }{
         Dec without E-terms (radians)
      }
      \sstsubsection{
         eq = double (Given)
      }{
         Besselian epoch of mean equator and equinox
      }
      \sstsubsection{
         rc = double $*$ (Returned)
      }{
         RA with E-terms included (radians)
      }
      \sstsubsection{
         dc = double $*$ (Returned)
      }{
         Dec with E-terms included (radians)
      }
   }
   \sstnotes{
      Most star positions from pre-1984 optical catalogues (or
      derived from astrometry using such stars) embody the
      E-terms.  If it is necessary to convert a formal mean
      place (for example a pulsar timing position) to one
      consistent with such a star catalogue, then the RA,Dec
      should be adjusted using this routine.
   }
   \sstdiytopic{
      See Also
   }{
      Explanatory Supplement to the Astronomical Ephemeris,
      section 2D, page 48.
   }
}
\sstroutine{
   palAirmas
}{
   Air mass at given zenith distance
}{
   \sstdescription{
      Calculates the airmass at the observed zenith distance.
   }
   \sstinvocation{
      double palAirmas( double zd );
   }
   \sstarguments{
      \sstsubsection{
         zd = double (Given)
      }{
         Observed zenith distance (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The \texttt{"} observed\texttt{"}  zenith distance referred to above means \texttt{"} as
           affected by refraction\texttt{"} .

         \sstitem
          Uses Hardie\texttt{'} s (1962) polynomial fit to Bemporad\texttt{'} s data for
           the relative air mass, X, in units of thickness at the zenith
           as tabulated by Schoenberg (1929). This is adequate for all
           normal needs as it is accurate to better than 0.1\% up to X =
           6.8 and better than 1\% up to X = 10. Bemporad\texttt{'} s tabulated
           values are unlikely to be trustworthy to such accuracy
           because of variations in density, pressure and other
           conditions in the atmosphere from those assumed in his work.

         \sstitem
          The sign of the ZD is ignored.

         \sstitem
          At zenith distances greater than about ZD = 87 degrees the
           air mass is held constant to avoid arithmetic overflows.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          Hardie, R.H., 1962, in \texttt{"} Astronomical Techniques\texttt{"}
             ed. W.A. Hiltner, University of Chicago Press, p180.

         \sstitem
          Schoenberg, E., 1929, Hdb. d. Ap.,
             Berlin, Julius Springer, 2, 268.
      }
   }
}
\sstroutine{
   palAltaz
}{
   Positions, velocities and accelerations for an altazimuth telescope mount
}{
   \sstdescription{
      Positions, velocities and accelerations for an altazimuth
      telescope mount.
   }
   \sstinvocation{
      palAltaz ( double ha, double dec, double phi,
                 double $*$az, double $*$azd, double $*$azdd,
                 double $*$el, double $*$eld, double $*$eldd,
                 double $*$pa, double $*$pad, double $*$padd );
   }
   \sstarguments{
      \sstsubsection{
         ha = double (Given)
      }{
         Hour angle (radians)
      }
      \sstsubsection{
         dec = double (Given)
      }{
         Declination (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observatory latitude (radians)
      }
      \sstsubsection{
         az = double $*$ (Returned)
      }{
         Azimuth (radians)
      }
      \sstsubsection{
         azd = double $*$ (Returned)
      }{
         Azimuth velocity (radians per radian of HA)
      }
      \sstsubsection{
         azdd = double $*$ (Returned)
      }{
         Azimuth acceleration (radians per radian of HA squared)
      }
      \sstsubsection{
         el = double $*$ (Returned)
      }{
         Elevation (radians)
      }
      \sstsubsection{
         eld = double $*$ (Returned)
      }{
         Elevation velocity (radians per radian of HA)
      }
      \sstsubsection{
         eldd = double $*$ (Returned)
      }{
         Elevation acceleration (radians per radian of HA squared)
      }
      \sstsubsection{
         pa = double $*$ (Returned)
      }{
         Parallactic angle (radians)
      }
      \sstsubsection{
         pad = double $*$ (Returned)
      }{
         Parallactic angle velocity (radians per radian of HA)
      }
      \sstsubsection{
         padd = double $*$ (Returned)
      }{
         Parallactic angle acceleration (radians per radian of HA squared)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Natural units are used throughout.  HA, DEC, PHI, AZ, EL
           and ZD are in radians.  The velocities and accelerations
           assume constant declination and constant rate of change of
           hour angle (as for tracking a star);  the units of AZD, ELD
           and PAD are radians per radian of HA, while the units of AZDD,
           ELDD and PADD are radians per radian of HA squared.  To
           convert into practical degree- and second-based units:

      }
          angles $*$ 360/2pi -$>$ degrees
          velocities $*$ (2pi/86400)$*$(360/2pi) -$>$ degree/sec
          accelerations $*$ ((2pi/86400)$*$$*$2)$*$(360/2pi) -$>$ degree/sec/sec

        Note that the seconds here are sidereal rather than SI.  One
        sidereal second is about 0.99727 SI seconds.

        The velocity and acceleration factors assume the sidereal
        tracking case.  Their respective numerical values are (exactly)
        1/240 and (approximately) 1/3300236.9.

      \sstitemlist{

         \sstitem
          Azimuth is returned in the range 0-2pi;  north is zero,
           and east is $+$pi/2.  Elevation and parallactic angle are
           returned in the range $+$/-pi.  Parallactic angle is $+$ve for
           a star west of the meridian and is the angle NP-star-zenith.

         \sstitem
          The latitude is geodetic as opposed to geocentric.  The
           hour angle and declination are topocentric.  Refraction and
           deficiencies in the telescope mounting are ignored.  The
           purpose of the routine is to give the general form of the
           quantities.  The details of a real telescope could profoundly
           change the results, especially close to the zenith.

         \sstitem
          No range checking of arguments is carried out.

         \sstitem
          In applications which involve many such calculations, rather
           than calling the present routine it will be more efficient to
           use inline code, having previously computed fixed terms such
           as sine and cosine of latitude, and (for tracking a star)
           sine and cosine of declination.
      }
   }
}
\sstroutine{
   palAmp
}{
   Convert star RA,Dec from geocentric apparaent to mean place
}{
   \sstdescription{
      Convert star RA,Dec from geocentric apparent to mean place. The
      mean coordinate system is close to ICRS. See palAmpqk for details.
   }
   \sstinvocation{
      void palAmp ( double ra, double da, double date, double eq,
                    double $*$rm, double $*$dm );
   }
   \sstarguments{
      \sstsubsection{
         ra = double (Given)
      }{
         Apparent RA (radians)
      }
      \sstsubsection{
         dec = double (Given)
      }{
         Apparent Dec (radians)
      }
      \sstsubsection{
         date = double (Given)
      }{
         TDB for apparent place (JD-2400000.5)
      }
      \sstsubsection{
         eq = double (Given)
      }{
         Equinox: Julian epoch of mean place.
      }
      \sstsubsection{
         rm = double $*$ (Returned)
      }{
         Mean RA (radians)
      }
      \sstsubsection{
         dm = double $*$ (Returned)
      }{
         Mean Dec (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          See palMappa and palAmpqk for details.
      }
   }
}
\sstroutine{
   palAmpqk
}{
   Convert star RA,Dec from geocentric apparent to mean place
}{
   \sstdescription{
      Convert star RA,Dec from geocentric apparent to mean place. The \texttt{"} mean\texttt{"}
      coordinate system is in fact close to ICRS. Use of this function
      is appropriate when efficiency is important and where many star
      positions are all to be transformed for one epoch and equinox.  The
      star-independent parameters can be obtained by calling the palMappa
      function.
   }
   \sstinvocation{
      void palAmpqk ( double ra, double da, double amprms[21],
                      double $*$rm, double $*$dm )
   }
   \sstarguments{
      \sstsubsection{
         ra = double (Given)
      }{
         Apparent RA (radians).
      }
      \sstsubsection{
         da = double (Given)
      }{
         Apparent Dec (radians).
      }
      \sstsubsection{
         amprms = double[21] (Given)
      }{
         Star-independent mean-to-apparent parameters (see palMappa):
         (0)      time interval for proper motion (Julian years)
         (1-3)    barycentric position of the Earth (AU)
         (4-6)    heliocentric direction of the Earth (unit vector)
         (7)      (grav rad Sun)$*$2/(Sun-Earth distance)
         (8-10)   abv: barycentric Earth velocity in units of c
         (11)     sqrt(1-v$*$v) where v=modulus(abv)
         (12-20)  precession/nutation (3,3) matrix
      }
      \sstsubsection{
         rm = double (Returned)
      }{
         Mean RA (radians).
      }
      \sstsubsection{
         dm = double (Returned)
      }{
         Mean Dec (radians).
      }
   }
   \sstdiytopic{
      Note
   }{
      Iterative techniques are used for the aberration and
      light deflection corrections so that the routines
      palAmp (or palAmpqk) and palMap (or palMapqk) are
      accurate inverses;  even at the edge of the Sun\texttt{'} s disc
      the discrepancy is only about 1 nanoarcsecond.
   }
}
\sstroutine{
   palAop
}{
   Apparent to observed place
}{
   \sstdescription{
      Apparent to observed place for sources distant from the solar system.
   }
   \sstinvocation{
      void palAop ( double rap, double dap, double date, double dut,
                    double elongm, double phim, double hm, double xp,
                    double yp, double tdk, double pmb, double rh,
                    double wl, double tlr,
                    double $*$aob, double $*$zob, double $*$hob,
                    double $*$dob, double $*$rob );
   }
   \sstarguments{
      \sstsubsection{
         rap = double (Given)
      }{
         Geocentric apparent right ascension
      }
      \sstsubsection{
         dap = double (Given)
      }{
         Geocentirc apparent declination
      }
      \sstsubsection{
         date = double (Given)
      }{
         UTC date/time (Modified Julian Date, JD-2400000.5)
      }
      \sstsubsection{
         dut = double (Given)
      }{
         delta UT: UT1-UTC (UTC seconds)
      }
      \sstsubsection{
         elongm = double (Given)
      }{
         Mean longitude of the observer (radians, east $+$ve)
      }
      \sstsubsection{
         phim = double (Given)
      }{
         Mean geodetic latitude of the observer (radians)
      }
      \sstsubsection{
         hm = double (Given)
      }{
         Observer\texttt{'} s height above sea level (metres)
      }
      \sstsubsection{
         xp = double (Given)
      }{
         Polar motion x-coordinates (radians)
      }
      \sstsubsection{
         yp = double (Given)
      }{
         Polar motion y-coordinates (radians)
      }
      \sstsubsection{
         tdk = double (Given)
      }{
         Local ambient temperature (K; std=273.15)
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         Local atmospheric pressure (mb; std=1013.25)
      }
      \sstsubsection{
         rh = double (Given)
      }{
         Local relative humidity (in the range 0.0-1.0)
      }
      \sstsubsection{
         wl = double (Given)
      }{
         Effective wavelength (micron, e.g. 0.55)
      }
      \sstsubsection{
         tlr = double (Given)
      }{
         Tropospheric laps rate (K/metre, e.g. 0.0065)
      }
      \sstsubsection{
         aob = double $*$ (Returned)
      }{
         Observed azimuth (radians: N=0; E=90)
      }
      \sstsubsection{
         zob = double $*$ (Returned)
      }{
         Observed zenith distance (radians)
      }
      \sstsubsection{
         hob = double $*$ (Returned)
      }{
         Observed Hour Angle (radians)
      }
      \sstsubsection{
         dob = double $*$ (Returned)
      }{
         Observed Declination (radians)
      }
      \sstsubsection{
         rob = double $*$ (Returned)
      }{
         Observed Right Ascension (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          This routine returns zenith distance rather than elevation
           in order to reflect the fact that no allowance is made for
           depression of the horizon.

         \sstitem
          The accuracy of the result is limited by the corrections for
           refraction.  Providing the meteorological parameters are
           known accurately and there are no gross local effects, the
           predicted apparent RA,Dec should be within about 0.1 arcsec
           for a zenith distance of less than 70 degrees.  Even at a
           topocentric zenith distance of 90 degrees, the accuracy in
           elevation should be better than 1 arcmin;  useful results
           are available for a further 3 degrees, beyond which the
           palRefro routine returns a fixed value of the refraction.
           The complementary routines palAop (or palAopqk) and palOap
           (or palOapqk) are self-consistent to better than 1 micro-
           arcsecond all over the celestial sphere.

         \sstitem
          It is advisable to take great care with units, as even
           unlikely values of the input parameters are accepted and
           processed in accordance with the models used.

         \sstitem
          \texttt{"} Apparent\texttt{"}  place means the geocentric apparent right ascension
           and declination, which is obtained from a catalogue mean place
           by allowing for space motion, parallax, precession, nutation,
           annual aberration, and the Sun\texttt{'} s gravitational lens effect.  For
           star positions in the FK5 system (i.e. J2000), these effects can
           be applied by means of the palMap etc routines.  Starting from
           other mean place systems, additional transformations will be
           needed;  for example, FK4 (i.e. B1950) mean places would first
           have to be converted to FK5, which can be done with the
           palFk425 etc routines.

         \sstitem
          \texttt{"} Observed\texttt{"}  Az,El means the position that would be seen by a
           perfect theodolite located at the observer.  This is obtained
           from the geocentric apparent RA,Dec by allowing for Earth
           orientation and diurnal aberration, rotating from equator
           to horizon coordinates, and then adjusting for refraction.
           The HA,Dec is obtained by rotating back into equatorial
           coordinates, using the geodetic latitude corrected for polar
           motion, and is the position that would be seen by a perfect
           equatorial located at the observer and with its polar axis
           aligned to the Earth\texttt{'} s axis of rotation (n.b. not to the
           refracted pole).  Finally, the RA is obtained by subtracting
           the HA from the local apparent ST.

         \sstitem
          To predict the required setting of a real telescope, the
           observed place produced by this routine would have to be
           adjusted for the tilt of the azimuth or polar axis of the
           mounting (with appropriate corrections for mount flexures),
           for non-perpendicularity between the mounting axes, for the
           position of the rotator axis and the pointing axis relative
           to it, for tube flexure, for gear and encoder errors, and
           finally for encoder zero points.  Some telescopes would, of
           course, exhibit other properties which would need to be
           accounted for at the appropriate point in the sequence.

         \sstitem
          This routine takes time to execute, due mainly to the
           rigorous integration used to evaluate the refraction.
           For processing multiple stars for one location and time,
           call palAoppa once followed by one call per star to palAopqk.
           Where a range of times within a limited period of a few hours
           is involved, and the highest precision is not required, call
           palAoppa once, followed by a call to palAoppat each time the
           time changes, followed by one call per star to palAopqk.

         \sstitem
          The DATE argument is UTC expressed as an MJD.  This is,
           strictly speaking, wrong, because of leap seconds.  However,
           as long as the delta UT and the UTC are consistent there
           are no difficulties, except during a leap second.  In this
           case, the start of the 61st second of the final minute should
           begin a new MJD day and the old pre-leap delta UT should
           continue to be used.  As the 61st second completes, the MJD
           should revert to the start of the day as, simultaneously,
           the delta UTC changes by one second to its post-leap new value.

         \sstitem
          The delta UT (UT1-UTC) is tabulated in IERS circulars and
           elsewhere.  It increases by exactly one second at the end of
           each UTC leap second, introduced in order to keep delta UT
           within $+$/- 0.9 seconds.

         \sstitem
          IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION.
           The longitude required by the present routine is east-positive,
           in accordance with geographical convention (and right-handed).
           In particular, note that the longitudes returned by the
           palObs routine are west-positive, following astronomical
           usage, and must be reversed in sign before use in the present
           routine.

         \sstitem
          The polar coordinates XP,YP can be obtained from IERS
           circulars and equivalent publications.  The maximum amplitude
           is about 0.3 arcseconds.  If XP,YP values are unavailable,
           use XP=YP=0.0.  See page B60 of the 1988 Astronomical Almanac
           for a definition of the two angles.

         \sstitem
          The height above sea level of the observing station, HM,
           can be obtained from the Astronomical Almanac (Section J
           in the 1988 edition), or via the routine palObs.  If P,
           the pressure in millibars, is available, an adequate
           estimate of HM can be obtained from the expression

      }
              HM $\sim$ -29.3$*$TSL$*$LOG(P/1013.25).

        where TSL is the approximate sea-level air temperature in K
        (see Astrophysical Quantities, C.W.Allen, 3rd edition,
        section 52).  Similarly, if the pressure P is not known,
        it can be estimated from the height of the observing
        station, HM, as follows:

              P $\sim$ 1013.25$*$EXP(-HM/(29.3$*$TSL)).

        Note, however, that the refraction is nearly proportional to the
        pressure and that an accurate P value is important for precise
        work.

      \sstitemlist{

         \sstitem
          The azimuths etc produced by the present routine are with
           respect to the celestial pole.  Corrections to the terrestrial
           pole can be computed using palPolmo.
      }
   }
}
\sstroutine{
   palAoppa
}{
   Precompute apparent to observed place parameters
}{
   \sstdescription{
      Precompute apparent to observed place parameters required by palAopqk
      and palOapqk.
   }
   \sstinvocation{
      void palAoppa ( double date, double dut, double elongm, double phim,
                      double hm, double xp, double yp, double tdk, double pmb,
                      double rh, double wl, double tlr, double aoprms[14] );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         UTC date/time (modified Julian Date, JD-2400000.5)
      }
      \sstsubsection{
         dut = double (Given)
      }{
         delta UT:  UT1-UTC (UTC seconds)
      }
      \sstsubsection{
         elongm = double (Given)
      }{
         mean longitude of the observer (radians, east $+$ve)
      }
      \sstsubsection{
         phim = double (Given)
      }{
         mean geodetic latitude of the observer (radians)
      }
      \sstsubsection{
         hm = double (Given)
      }{
         observer\texttt{'} s height above sea level (metres)
      }
      \sstsubsection{
         xp = double (Given)
      }{
         polar motion x-coordinate (radians)
      }
      \sstsubsection{
         yp = double (Given)
      }{
         polar motion y-coordinate (radians)
      }
      \sstsubsection{
         tdk = double (Given)
      }{
         local ambient temperature (K; std=273.15)
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         local atmospheric pressure (mb; std=1013.25)
      }
      \sstsubsection{
         rh = double (Given)
      }{
         local relative humidity (in the range 0.0-1.0)
      }
      \sstsubsection{
         wl = double (Given)
      }{
         effective wavelength (micron, e.g. 0.55)
      }
      \sstsubsection{
         tlr = double (Given)
      }{
         tropospheric lapse rate (K/metre, e.g. 0.0065)
      }
      \sstsubsection{
         aoprms = double [14] (Returned)
      }{
         Star-independent apparent-to-observed parameters

          (0)      geodetic latitude (radians)
          (1,2)    sine and cosine of geodetic latitude
          (3)      magnitude of diurnal aberration vector
          (4)      height (hm)
          (5)      ambient temperature (tdk)
          (6)      pressure (pmb)
          (7)      relative humidity (rh)
          (8)      wavelength (wl)
          (9)     lapse rate (tlr)
          (10,11)  refraction constants A and B (radians)
          (12)     longitude $+$ eqn of equinoxes $+$ sidereal DUT (radians)
          (13)     local apparent sidereal time (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          It is advisable to take great care with units, as even
           unlikely values of the input parameters are accepted and
           processed in accordance with the models used.

         \sstitem
          The DATE argument is UTC expressed as an MJD.  This is,
           strictly speaking, improper, because of leap seconds.  However,
           as long as the delta UT and the UTC are consistent there
           are no difficulties, except during a leap second.  In this
           case, the start of the 61st second of the final minute should
           begin a new MJD day and the old pre-leap delta UT should
           continue to be used.  As the 61st second completes, the MJD
           should revert to the start of the day as, simultaneously,
           the delta UTC changes by one second to its post-leap new value.

         \sstitem
          The delta UT (UT1-UTC) is tabulated in IERS circulars and
           elsewhere.  It increases by exactly one second at the end of
           each UTC leap second, introduced in order to keep delta UT
           within $+$/- 0.9 seconds.

         \sstitem
          IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION.
           The longitude required by the present routine is east-positive,
           in accordance with geographical convention (and right-handed).
           In particular, note that the longitudes returned by the
           palObs routine are west-positive, following astronomical
           usage, and must be reversed in sign before use in the present
           routine.

         \sstitem
          The polar coordinates XP,YP can be obtained from IERS
           circulars and equivalent publications.  The maximum amplitude
           is about 0.3 arcseconds.  If XP,YP values are unavailable,
           use XP=YP=0.0.  See page B60 of the 1988 Astronomical Almanac
           for a definition of the two angles.

         \sstitem
          The height above sea level of the observing station, HM,
           can be obtained from the Astronomical Almanac (Section J
           in the 1988 edition), or via the routine palObs.  If P,
           the pressure in millibars, is available, an adequate
           estimate of HM can be obtained from the expression

      }
              HM $\sim$ -29.3$*$TSL$*$log(P/1013.25).

        where TSL is the approximate sea-level air temperature in K
        (see Astrophysical Quantities, C.W.Allen, 3rd edition,
        section 52).  Similarly, if the pressure P is not known,
        it can be estimated from the height of the observing
        station, HM, as follows:

              P $\sim$ 1013.25$*$exp(-HM/(29.3$*$TSL)).

        Note, however, that the refraction is nearly proportional to the
        pressure and that an accurate P value is important for precise
        work.

      \sstitemlist{

         \sstitem
          Repeated, computationally-expensive, calls to palAoppa for
           times that are very close together can be avoided by calling
           palAoppa just once and then using palAoppat for the subsequent
           times.  Fresh calls to palAoppa will be needed only when
           changes in the precession have grown to unacceptable levels or
           when anything affecting the refraction has changed.
      }
   }
}
\sstroutine{
   palAoppat
}{
   Recompute sidereal time to support apparent to observed place
}{
   \sstdescription{
      This routine recomputes the sidereal time in the apparent to
      observed place star-independent parameter block.
   }
   \sstinvocation{
      void palAoppat( double date, double aoprms[14] );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         UTC date/time (modified Julian Date, JD-2400000.5)
         (see palAoppa description for comments on leap seconds)
      }
      \sstsubsection{
         aoprms = double[14] (Given \& Returned)
      }{
         Star-independent apparent-to-observed parameters. Updated
         by this routine. Requires element 12 to be the longitude $+$
         eqn of equinoxes $+$ sidereal DUT and fills in element 13
         with the local apparent sidereal time (in radians).
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          See palAoppa for more information.

         \sstitem
          The star-independent parameters are not treated as an opaque
           struct in order to retain compatibility with SLA.
      }
   }
}
\sstroutine{
   palAopqk
}{
   Quick apparent to observed place
}{
   \sstdescription{
      Quick apparent to observed place.
   }
   \sstinvocation{
      void palAopqk ( double rap, double dap, const double aoprms[14],
                      double $*$aob, double $*$zob, double $*$hob,
                      double $*$dob, double $*$rob );
   }
   \sstarguments{
      \sstsubsection{
         rap = double (Given)
      }{
         Geocentric apparent right ascension
      }
      \sstsubsection{
         dap = double (Given)
      }{
         Geocentric apparent declination
      }
      \sstsubsection{
         aoprms = const double [14] (Given)
      }{
         Star-independent apparent-to-observed parameters.

          [0]      geodetic latitude (radians)
          [1,2]    sine and cosine of geodetic latitude
          [3]      magnitude of diurnal aberration vector
          [4]      height (HM)
          [5]      ambient temperature (T)
          [6]      pressure (P)
          [7]      relative humidity (RH)
          [8]      wavelength (WL)
          [9]      lapse rate (TLR)
          [10,11]  refraction constants A and B (radians)
          [12]     longitude $+$ eqn of equinoxes $+$ sidereal DUT (radians)
          [13]     local apparent sidereal time (radians)
      }
      \sstsubsection{
         aob = double $*$ (Returned)
      }{
         Observed azimuth (radians: N=0,E=90)
      }
      \sstsubsection{
         zob = double $*$ (Returned)
      }{
         Observed zenith distance (radians)
      }
      \sstsubsection{
         hob = double $*$ (Returned)
      }{
         Observed Hour Angle (radians)
      }
      \sstsubsection{
         dob = double $*$ (Returned)
      }{
         Observed Declination (radians)
      }
      \sstsubsection{
         rob = double $*$ (Returned)
      }{
         Observed Right Ascension (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          This routine returns zenith distance rather than elevation
           in order to reflect the fact that no allowance is made for
           depression of the horizon.

         \sstitem
          The accuracy of the result is limited by the corrections for
           refraction.  Providing the meteorological parameters are
           known accurately and there are no gross local effects, the
           observed RA,Dec predicted by this routine should be within
           about 0.1 arcsec for a zenith distance of less than 70 degrees.
           Even at a topocentric zenith distance of 90 degrees, the
           accuracy in elevation should be better than 1 arcmin;  useful
           results are available for a further 3 degrees, beyond which
           the palRefro routine returns a fixed value of the refraction.
           The complementary routines palAop (or palAopqk) and palOap
           (or palOapqk) are self-consistent to better than 1 micro-
           arcsecond all over the celestial sphere.

         \sstitem
          It is advisable to take great care with units, as even
           unlikely values of the input parameters are accepted and
           processed in accordance with the models used.

         \sstitem
          \texttt{"} Apparent\texttt{"}  place means the geocentric apparent right ascension
           and declination, which is obtained from a catalogue mean place
           by allowing for space motion, parallax, precession, nutation,
           annual aberration, and the Sun\texttt{'} s gravitational lens effect.  For
           star positions in the FK5 system (i.e. J2000), these effects can
           be applied by means of the palMap etc routines.  Starting from
           other mean place systems, additional transformations will be
           needed;  for example, FK4 (i.e. B1950) mean places would first
           have to be converted to FK5, which can be done with the
           palFk425 etc routines.

         \sstitem
          \texttt{"} Observed\texttt{"}  Az,El means the position that would be seen by a
           perfect theodolite located at the observer.  This is obtained
           from the geocentric apparent RA,Dec by allowing for Earth
           orientation and diurnal aberration, rotating from equator
           to horizon coordinates, and then adjusting for refraction.
           The HA,Dec is obtained by rotating back into equatorial
           coordinates, using the geodetic latitude corrected for polar
           motion, and is the position that would be seen by a perfect
           equatorial located at the observer and with its polar axis
           aligned to the Earth\texttt{'} s axis of rotation (n.b. not to the
           refracted pole).  Finally, the RA is obtained by subtracting
           the HA from the local apparent ST.

         \sstitem
          To predict the required setting of a real telescope, the
           observed place produced by this routine would have to be
           adjusted for the tilt of the azimuth or polar axis of the
           mounting (with appropriate corrections for mount flexures),
           for non-perpendicularity between the mounting axes, for the
           position of the rotator axis and the pointing axis relative
           to it, for tube flexure, for gear and encoder errors, and
           finally for encoder zero points.  Some telescopes would, of
           course, exhibit other properties which would need to be
           accounted for at the appropriate point in the sequence.

         \sstitem
          The star-independent apparent-to-observed-place parameters
           in AOPRMS may be computed by means of the palAoppa routine.
           If nothing has changed significantly except the time, the
           palAoppat routine may be used to perform the requisite
           partial recomputation of AOPRMS.

         \sstitem
          At zenith distances beyond about 76 degrees, the need for
           special care with the corrections for refraction causes a
           marked increase in execution time.  Moreover, the effect
           gets worse with increasing zenith distance.  Adroit
           programming in the calling application may allow the
           problem to be reduced.  Prepare an alternative AOPRMS array,
           computed for zero air-pressure;  this will disable the
           refraction corrections and cause rapid execution.  Using
           this AOPRMS array, a preliminary call to the present routine
           will, depending on the application, produce a rough position
           which may be enough to establish whether the full, slow
           calculation (using the real AOPRMS array) is worthwhile.
           For example, there would be no need for the full calculation
           if the preliminary call had already established that the
           source was well below the elevation limits for a particular
           telescope.

         \sstitem
          The azimuths etc produced by the present routine are with
           respect to the celestial pole.  Corrections to the terrestrial
           pole can be computed using palPolmo.
      }
   }
}
\sstroutine{
   palAtmdsp
}{
   Apply atmospheric-dispersion adjustments to refraction coefficients
}{
   \sstdescription{
      Apply atmospheric-dispersion adjustments to refraction coefficients.
   }
   \sstinvocation{
      void palAtmdsp( double tdk, double pmb, double rh, double wl1,
                      double a1, double b1, double wl2, double $*$a2, double $*$b2 );
   }
   \sstarguments{
      \sstsubsection{
         tdk = double (Given)
      }{
         Ambient temperature, K
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         Ambient pressure, millibars
      }
      \sstsubsection{
         rh = double (Given)
      }{
         Ambient relative humidity, 0-1
      }
      \sstsubsection{
         wl1 = double (Given)
      }{
         Reference wavelength, micrometre (0.4 recommended)
      }
      \sstsubsection{
         a1 = double (Given)
      }{
         Refraction coefficient A for wavelength wl1 (radians)
      }
      \sstsubsection{
         b1 = double (Given)
      }{
         Refraction coefficient B for wavelength wl1 (radians)
      }
      \sstsubsection{
         wl2 = double (Given)
      }{
         Wavelength for which adjusted A,B required
      }
      \sstsubsection{
         a2 = double $*$ (Returned)
      }{
         Refraction coefficient A for wavelength WL2 (radians)
      }
      \sstsubsection{
         b2 = double $*$ (Returned)
      }{
         Refraction coefficient B for wavelength WL2 (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          To use this routine, first call palRefco specifying WL1 as the
         wavelength.  This yields refraction coefficients A1,B1, correct
         for that wavelength.  Subsequently, calls to palAtmdsp specifying
         different wavelengths will produce new, slightly adjusted
         refraction coefficients which apply to the specified wavelength.

         \sstitem
          Most of the atmospheric dispersion happens between 0.7 micrometre
         and the UV atmospheric cutoff, and the effect increases strongly
         towards the UV end.  For this reason a blue reference wavelength
         is recommended, for example 0.4 micrometres.

         \sstitem
          The accuracy, for this set of conditions:

      }
         height above sea level    2000 m
                       latitude    29 deg
                       pressure    793 mb
                    temperature    17 degC
                       humidity    50\%
                     lapse rate    0.0065 degC/m
           reference wavelength    0.4 micrometre
                 star elevation    15 deg

      is about 2.5 mas RMS between 0.3 and 1.0 micrometres, and stays
      within 4 mas for the whole range longward of 0.3 micrometres
      (compared with a total dispersion from 0.3 to 20.0 micrometres
      of about 11 arcsec).  These errors are typical for ordinary
      conditions and the given elevation;  in extreme conditions values
      a few times this size may occur, while at higher elevations the
      errors become much smaller.

      \sstitemlist{

         \sstitem
          If either wavelength exceeds 100 micrometres, the radio case
         is assumed and the returned refraction coefficients are the
         same as the given ones.  Note that radio refraction coefficients
         cannot be turned into optical values using this routine, nor
         vice versa.

         \sstitem
          The algorithm consists of calculation of the refractivity of the
         air at the observer for the two wavelengths, using the methods
         of the palRefro routine, and then scaling of the two refraction
         coefficients according to classical refraction theory.  This
         amounts to scaling the A coefficient in proportion to (n-1) and
         the B coefficient almost in the same ratio (see R.M.Green,
         \texttt{"} Spherical Astronomy\texttt{"} , Cambridge University Press, 1985).
      }
   }
}
\sstroutine{
   palCaldj
}{
   Gregorian Calendar to Modified Julian Date
}{
   \sstdescription{
      Modified Julian Date to Gregorian Calendar with special
      behaviour for 2-digit years relating to 1950 to 2049.
   }
   \sstinvocation{
      void palCaldj ( int iy, int im, int id, double $*$djm, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         iy = int (Given)
      }{
         Year in the Gregorian calendar
      }
      \sstsubsection{
         im = int (Given)
      }{
         Month in the Gergorian calendar
      }
      \sstsubsection{
         id = int (Given)
      }{
         Day in the Gregorian calendar
      }
      \sstsubsection{
         djm = double $*$ (Returned)
      }{
         Modified Julian Date (JD-2400000.5) for 0 hrs
      }
      \sstsubsection{
         j = status (Returned)
      }{
         0 = OK. See eraCal2jd for other values.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraCal2jd

         \sstitem
          Unlike eraCal2jd this routine treats the years 0-100 as
           referring to the end of the 20th Century and beginning of
           the 21st Century. If this behaviour is not acceptable
           use the SOFA/ERFA routine directly or palCldj.
           Acceptable years are 00-49, interpreted as 2000-2049,
                                50-99,     \texttt{"}        \texttt{"}   1950-1999,
                                all others, interpreted literally.

         \sstitem
          Unlike SLA this routine will work with negative years.
      }
   }
}
\sstroutine{
   palDafin
}{
   Sexagesimal character string to angle
}{
   \sstdescription{
      Extracts an angle from a sexagesimal string with degrees, arcmin,
      arcsec fields using space or comma delimiters.
   }
   \sstinvocation{
      void palDafin ( const char $*$string, int $*$ipos, double $*$a, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         string = const char $*$ (Given)
      }{
         String containing deg, arcmin, arcsec fields
      }
      \sstsubsection{
         ipos = int $*$ (Given \& Returned)
      }{
         Position to start decoding \texttt{"} string\texttt{"} . First character
         is position 1 for compatibility with SLA. After
         calling this routine \texttt{"} iptr\texttt{"}  will be positioned after
         the sexagesimal string.
      }
      \sstsubsection{
         a = double $*$ (Returned)
      }{
         Angle in radians.
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         status:  0 = OK
                 $+$1 = default, A unchanged
         \sstitemlist{

            \sstitem
                     1 = bad degrees      )

            \sstitem
                     2 = bad arcminutes   )  (note 3)

            \sstitem
                     3 = bad arcseconds   )
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The first three \texttt{"} fields\texttt{"}  in STRING are degrees, arcminutes,
           arcseconds, separated by spaces or commas.  The degrees field
           may be signed, but not the others.  The decoding is carried
           out by the palDfltin routine and is free-format.

         \sstitem
          Successive fields may be absent, defaulting to zero.  For
           zero status, the only combinations allowed are degrees alone,
           degrees and arcminutes, and all three fields present.  If all
           three fields are omitted, a status of $+$1 is returned and A is
           unchanged.  In all other cases A is changed.

         \sstitem
          Range checking:

      }
            The degrees field is not range checked.  However, it is
            expected to be integral unless the other two fields are absent.

            The arcminutes field is expected to be 0-59, and integral if
            the arcseconds field is present.  If the arcseconds field
            is absent, the arcminutes is expected to be 0-59.9999...

            The arcseconds field is expected to be 0-59.9999...

      \sstitemlist{

         \sstitem
          Decoding continues even when a check has failed.  Under these
           circumstances the field takes the supplied value, defaulting
           to zero, and the result A is computed and returned.

         \sstitem
          Further fields after the three expected ones are not treated
           as an error.  The pointer IPOS is left in the correct state
           for further decoding with the present routine or with palDfltin
           etc. See the example, above.

         \sstitem
          If STRING contains hours, minutes, seconds instead of degrees
           etc, or if the required units are turns (or days) instead of
           radians, the result A should be multiplied as follows:

      }
            for        to obtain    multiply
            STRING     A in         A by

            d \texttt{'}  \texttt{"}       radians      1       =  1.0
            d \texttt{'}  \texttt{"}       turns        1/2pi   =  0.1591549430918953358
            h m s      radians      15      =  15.0
            h m s      days         15/2pi  =  2.3873241463784300365
   }
   \sstdiytopic{
      Example
   }{
      argument    before                           after

      STRING      \texttt{'} -57 17 44.806  12 34 56.7\texttt{'}       unchanged
      IPTR        1                                16 (points to 12...)
      A           ?                                -1.00000D0
      J           ?                                0
   }
}
\sstroutine{
   palDcmpf
}{
   Decompose an [x,y] linear fit into its constituent parameters:
   zero points, scales, nonperpendicularity and orientation
}{
   \sstdescription{
      The model relates two sets of [x,y] coordinates as follows.
      Naming the elements of coeffs:
      \newline
      \newline
      \hspace*{1.5 em}
         coeffs[0] = A
      \newline
      \hspace*{1.5 em}
         coeffs[1] = B
      \newline
      \hspace*{1.5 em}
         coeffs[2] = C
      \newline
      \hspace*{1.5 em}
         coeffs[3] = D
      \newline
      \hspace*{1.5 em}
         coeffs[4] = E
      \newline
      \hspace*{1.5 em}
         coeffs[5] = F
      \newline
      \newline
      the model transforms coordinates [x1,y1] into coordinates
      [x2,y2] as follows:
      \newline
      \newline
      \hspace*{1.5 em}
         x2 = A $+$ B $*$ x1 $+$ C $*$ y1
      \newline
      \hspace*{1.5 em}
         y2 = D $+$ E $*$ x1 $+$ F $*$ y1
      \newline
      \newline
      The transformation can be decomposed into four steps:
      \newline
      \newline
      1) Zero points:
      \newline
      \hspace*{1.5 em}
         x\texttt{'}  = xz $+$ x1
      \newline
      \hspace*{1.5 em}
         y\texttt{'}  = yz $+$ y1
      \newline
      \newline
      2) Scales:
      \newline
      \hspace*{1.5 em}
         x\texttt{'} \texttt{'}  = xs $*$ x\texttt{'}
      \newline
      \hspace*{1.5 em}
         y\texttt{'} \texttt{'}  = ys $*$ y\texttt{'}
      \newline
      \newline
      3) Nonperpendicularity:
      \newline
      \hspace*{1.5 em}
         x\texttt{'} \texttt{'} \texttt{'}  = cos(perp / 2) $*$ x\texttt{'} \texttt{'}  $+$ sin(perp / 2) $*$ y\texttt{'} \texttt{'}
      \newline
      \hspace*{1.5 em}
         y\texttt{'} \texttt{'} \texttt{'}  = sin(perp / 2) $*$ x\texttt{'} \texttt{'}  $+$ cos(perp / 2) $*$ y\texttt{'} \texttt{'}
      \newline
      \newline
      4) Orientation:
      \newline
      \hspace*{1.5 em}
         x2 =  cos(orient) $*$ x\texttt{'} \texttt{'} \texttt{'}  $+$ sin(orient) $*$ y\texttt{'} \texttt{'} \texttt{'}
      \newline
      \hspace*{1.5 em}
         y2 = -sin(orient) $*$ y\texttt{'} \texttt{'} \texttt{'}  $+$ cos(orient) $*$ y\texttt{'} \texttt{'} \texttt{'}
      \newline
      \newline
   }
   \sstinvocation{
      palDcmpf ( double coeffs[6], double $*$xz, double $*$yz, double $*$xs,
                 double $*$ys, double $*$perp, double $*$orient )
   }
   \sstarguments{
      \sstsubsection{
         coeffs = double[6] (Given)
      }{
         transformation coefficients (see note)
      }
      \sstsubsection{
         xz = double (Returned)
      }{
         x zero point
      }
      \sstsubsection{
         yz = double (Returned)
      }{
         y zero point
      }
      \sstsubsection{
         xs = double (Returned)
      }{
         x scale
      }
      \sstsubsection{
         ys = double (Returned)
      }{
         y scale
      }
      \sstsubsection{
         perp = double (Returned)
      }{
         nonperpendicularity (radians)
      }
      \sstsubsection{
         orient = double (Returned)
      }{
         orientation (radians)
      }
   }
   \sstdiytopic{
      See also
   }{
      palFitxy, palPxy, palInvf and palXy2xy
   }
}
\sstroutine{
   palDe2h
}{
   Equatorial to horizon coordinates: HA,Dec to Az,E
}{
   \sstdescription{
      Convert equatorial to horizon coordinates.
   }
   \sstinvocation{
      palDe2h( double ha, double dec, double phi, double $*$ az, double $*$ el );
   }
   \sstarguments{
      \sstsubsection{
         ha = double $*$ (Given)
      }{
         Hour angle (radians)
      }
      \sstsubsection{
         dec = double $*$ (Given)
      }{
         Declination (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observatory latitude (radians)
      }
      \sstsubsection{
         az = double $*$ (Returned)
      }{
         Azimuth (radians)
      }
      \sstsubsection{
         el = double $*$ (Returned)
      }{
         Elevation (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          All the arguments are angles in radians.

         \sstitem
          Azimuth is returned in the range 0-2pi;  north is zero,
           and east is $+$pi/2.  Elevation is returned in the range
           $+$/-pi/2.

         \sstitem
          The latitude must be geodetic.  In critical applications,
           corrections for polar motion should be applied.

         \sstitem
          In some applications it will be important to specify the
           correct type of hour angle and declination in order to
           produce the required type of azimuth and elevation.  In
           particular, it may be important to distinguish between
           elevation as affected by refraction, which would
           require the \texttt{"} observed\texttt{"}  HA,Dec, and the elevation
           in vacuo, which would require the \texttt{"} topocentric\texttt{"}  HA,Dec.
           If the effects of diurnal aberration can be neglected, the
           \texttt{"} apparent\texttt{"}  HA,Dec may be used instead of the topocentric
           HA,Dec.

         \sstitem
          No range checking of arguments is carried out.

         \sstitem
          In applications which involve many such calculations, rather
           than calling the present routine it will be more efficient to
           use inline code, having previously computed fixed terms such
           as sine and cosine of latitude, and (for tracking a star)
           sine and cosine of declination.
      }
   }
}
\sstroutine{
   palDeuler
}{
   Form a rotation matrix from the Euler angles
}{
   \sstdescription{
      A rotation is positive when the reference frame rotates
      anticlockwise as seen looking towards the origin from the
      positive region of the specified axis.

      The characters of ORDER define which axes the three successive
      rotations are about.  A typical value is \texttt{'} ZXZ\texttt{'} , indicating that
      RMAT is to become the direction cosine matrix corresponding to
      rotations of the reference frame through PHI radians about the
      old Z-axis, followed by THETA radians about the resulting X-axis,
      then PSI radians about the resulting Z-axis.

      The axis names can be any of the following, in any order or
      combination:  X, Y, Z, uppercase or lowercase, 1, 2, 3.  Normal
      axis labelling/numbering conventions apply;  the xyz (=123)
      triad is right-handed.  Thus, the \texttt{'} ZXZ\texttt{'}  example given above
      could be written \texttt{'} zxz\texttt{'}  or \texttt{'} 313\texttt{'}  (or even \texttt{'} ZxZ\texttt{'}  or \texttt{'} 3xZ\texttt{'} ).  ORDER
      is terminated by length or by the first unrecognized character.

      Fewer than three rotations are acceptable, in which case the later
      angle arguments are ignored.  If all rotations are zero, the
      identity matrix is produced.
   }
   \sstinvocation{
      void palDeuler ( const char $*$order, double phi, double theta, double psi,
                       double rmat[3][3] );
   }
   \sstarguments{
      \sstsubsection{
         order = const char[] (Given)
      }{
         Specifies about which axes the rotation occurs
      }
      \sstsubsection{
         phi = double (Given)
      }{
         1st rotation (radians)
      }
      \sstsubsection{
         theta = double (Given)
      }{
         2nd rotation (radians)
      }
      \sstsubsection{
         psi = double (Given)
      }{
         3rd rotation (radians)
      }
      \sstsubsection{
         rmat = double[3][3] (Given \& Returned)
      }{
         Rotation matrix
      }
   }
}
\sstroutine{
   palDfltin
}{
   Convert free-format input into double precision floating point
}{
   \sstdescription{
      Extracts a number from an input string starting at the specified
      index.
   }
   \sstinvocation{
      void palDfltin( const char $*$ string, int $*$nstrt,
                      double $*$dreslt, int $*$jflag );
   }
   \sstarguments{
      \sstsubsection{
         string = const char $*$ (Given)
      }{
         String containing number to be decoded.
      }
      \sstsubsection{
         nstrt = int $*$ (Given and Returned)
      }{
         Character number indicating where decoding should start.
         On output its value is updated to be the location of the
         possible next value. For compatibility with SLA the first
         character is index 1.
      }
      \sstsubsection{
         dreslt = double $*$ (Returned)
      }{
         Result. Not updated when jflag=1.
      }
      \sstsubsection{
         jflag = int $*$ (Returned)
      }{
         status: -1 = -OK, 0 = $+$OK, 1 = null, 2 = error
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses the strtod() system call to do the parsing. This may lead to
           subtle differences when compared to the SLA/F parsing.

         \sstitem
          All \texttt{"} D\texttt{"}  characters are converted to \texttt{"} E\texttt{"}  to handle fortran exponents.

         \sstitem
          Commas are recognized as a special case and are skipped if one happens
           to be the next character when updating nstrt. Additionally the output
           nstrt position will skip past any trailing space.

         \sstitem
          If no number can be found flag will be set to 1.

         \sstitem
          If the number overflows or underflows jflag will be set to 2. For overflow
           the returned result will have the value HUGE\_VAL, for underflow it
           will have the value 0.0.

         \sstitem
          For compatiblity with SLA/F -0 will be returned as \texttt{"} 0\texttt{"}  with jflag == -1.

         \sstitem
          Unlike slaDfltin a standalone \texttt{"} E\texttt{"}  will return status 1 (could not find
           a number) rather than 2 (bad number).
      }
   }
   \sstimplementationstatus{
      \sstitemlist{

         \sstitem
          The code is more robust if the C99 copysign() function is available.
         This can recognize the -0.0 values returned by strtod. If copysign() is
         missing we try to scan the string looking for minus signs.
      }
   }
}
\sstroutine{
   palDh2e
}{
   Horizon to equatorial coordinates: Az,El to HA,Dec
}{
   \sstdescription{
      Convert horizon to equatorial coordinates.
   }
   \sstinvocation{
      palDh2e( double az, double el, double phi, double $*$ ha, double $*$ dec );
   }
   \sstarguments{
      \sstsubsection{
         az = double (Given)
      }{
         Azimuth (radians)
      }
      \sstsubsection{
         el = double (Given)
      }{
         Elevation (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observatory latitude (radians)
      }
      \sstsubsection{
         ha = double $*$ (Returned)
      }{
         Hour angle (radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Declination (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          All the arguments are angles in radians.

         \sstitem
          The sign convention for azimuth is north zero, east $+$pi/2.

         \sstitem
          HA is returned in the range $+$/-pi.  Declination is returned
           in the range $+$/-pi/2.

         \sstitem
          The latitude is (in principle) geodetic.  In critical
           applications, corrections for polar motion should be applied.

         \sstitem
          In some applications it will be important to specify the
           correct type of elevation in order to produce the required
           type of HA,Dec.  In particular, it may be important to
           distinguish between the elevation as affected by refraction,
           which will yield the \texttt{"} observed\texttt{"}  HA,Dec, and the elevation
           in vacuo, which will yield the \texttt{"} topocentric\texttt{"}  HA,Dec.  If the
           effects of diurnal aberration can be neglected, the
           topocentric HA,Dec may be used as an approximation to the
           \texttt{"} apparent\texttt{"}  HA,Dec.

         \sstitem
          No range checking of arguments is done.

         \sstitem
          In applications which involve many such calculations, rather
           than calling the present routine it will be more efficient to
           use inline code, having previously computed fixed terms such
           as sine and cosine of latitude.
      }
   }
}
\sstroutine{
   palDjcal
}{
   Modified Julian Date to Gregorian Calendar
}{
   \sstdescription{
      Modified Julian Date to Gregorian Calendar, expressed
      in a form convenient for formatting messages (namely
      rounded to a specified precision, and with the fields
      stored in a single array)
   }
   \sstinvocation{
      void palDjcal ( int ndp, double djm, int iymdf[4], int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         ndp = int (Given)
      }{
         Number of decimal places of days in fraction.
      }
      \sstsubsection{
         djm = double (Given)
      }{
         Modified Julian Date (JD-2400000.5)
      }
      \sstsubsection{
         iymdf[4] = int[] (Returned)
      }{
         Year, month, day, fraction in Gregorian calendar.
      }
      \sstsubsection{
         j = status (Returned)
      }{
         0 = OK. See eraJd2cal for other values.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraJd2cal
      }
   }
}
\sstroutine{
   palDmat
}{
   Matrix inversion \& solution of simultaneous equations
}{
   \sstdescription{
      Matrix inversion \& solution of simultaneous equations
      For the set of n simultaneous equations in n unknowns:
           A.Y = X
      this routine calculates the inverse of A, the determinant
      of matrix A and the vector of N unknowns.
   }
   \sstinvocation{
      void palDmat( int n, double $*$a, double $*$y, double $*$d, int $*$jf,
                     int $*$iw );
   }
   \sstarguments{
      \sstsubsection{
         n = int (Given)
      }{
         Number of simultaneous equations and number of unknowns.
      }
      \sstsubsection{
         a = double[] (Given \& Returned)
      }{
         A non-singular NxN matrix (implemented as a contiguous block
         of memory). After calling this routine \texttt{"} a\texttt{"}  contains the
         inverse of the matrix.
      }
      \sstsubsection{
         y = double[] (Given \& Returned)
      }{
         On input the vector of N knowns. On exit this vector contains the
         N solutions.
      }
      \sstsubsection{
         d = double $*$ (Returned)
      }{
         The determinant.
      }
      \sstsubsection{
         jf = int $*$ (Returned)
      }{
         The singularity flag.  If the matrix is non-singular, jf=0
         is returned.  If the matrix is singular, jf=-1 \& d=0.0 are
         returned.  In the latter case, the contents of array \texttt{"} a\texttt{"}  on
         return are undefined.
      }
      \sstsubsection{
         iw = int[] (Given)
      }{
         Integer workspace of size N.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Implemented using Gaussian elimination with partial pivoting.

         \sstitem
          Optimized for speed rather than accuracy with errors 1 to 4
           times those of routines optimized for accuracy.
      }
   }
}
\sstroutine{
   palDs2tp
}{
   Spherical to tangent plane projection
}{
   \sstdescription{
      Projection of spherical coordinates onto tangent plane:
      \texttt{"} gnomonic\texttt{"}  projection - \texttt{"} standard coordinates\texttt{"}
   }
   \sstinvocation{
      palDs2tp( double ra, double dec, double raz, double decz,
                double $*$xi, double $*$eta, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         ra = double (Given)
      }{
         RA spherical coordinate of point to be projected (radians)
      }
      \sstsubsection{
         dec = double (Given)
      }{
         Dec spherical coordinate of point to be projected (radians)
      }
      \sstsubsection{
         raz = double (Given)
      }{
         RA spherical coordinate of tangent point (radians)
      }
      \sstsubsection{
         decz = double (Given)
      }{
         Dec spherical coordinate of tangent point (radians)
      }
      \sstsubsection{
         xi = double $*$ (Returned)
      }{
         First rectangular coordinate on tangent plane (radians)
      }
      \sstsubsection{
         eta = double $*$ (Returned)
      }{
         Second rectangular coordinate on tangent plane (radians)
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         status: 0 = OK, star on tangent plane
                 1 = error, star too far from axis
                 2 = error, antistar on tangent plane
                 3 = error, antistar too far from axis
      }
   }
}
\sstroutine{
   palDat
}{
   Return offset between UTC and TAI
}{
   \sstdescription{
      Increment to be applied to Coordinated Universal Time UTC to give
      International Atomic Time (TAI).
   }
   \sstinvocation{
      dat = palDat( double utc );
   }
   \sstarguments{
      \sstsubsection{
         utc = double (Given)
      }{
         UTC date as a modified JD (JD-2400000.5)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         dat = double
      }{
         TAI-UTC in seconds
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          This routine converts the MJD argument to calendar date before calling
           the SOFA/ERFA eraDat function.

         \sstitem
          This routine matches the slaDat interface which differs from the eraDat
           interface. Consider coding directly to the SOFA/ERFA interface.

         \sstitem
          See eraDat for a description of error conditions when calling this function
           with a time outside of the UTC range.

         \sstitem
          The status argument from eraDat is ignored. This is reasonable since the
           error codes are mainly related to incorrect calendar dates when calculating
           the JD internally.
      }
   }
}
\sstroutine{
   palDmoon
}{
   Approximate geocentric position and velocity of the Moon
}{
   \sstdescription{
      Calculate the approximate geocentric position of the Moon
      using a full implementation of the algorithm published by
      Meeus (l\texttt{'} Astronomie, June 1984, p348).
   }
   \sstinvocation{
      void palDmoon( double date, double pv[6] );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TDB as a Modified Julian Date (JD-2400000.5)
      }
      \sstsubsection{
         pv = double [6] (Returned)
      }{
         Moon x,y,z,xdot,ydot,zdot, mean equator and
         equinox of date (AU, AU/s)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in
           latitude and 0.2 arcsec in HP (equivalent to about 20 km in
           distance).  Comparison with JPL DE200 over the interval
           1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in
           longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km
           and 81 mm/s in distance.  The maximum errors over the same
           interval are 18 arcsec and 0.50 arcsec/hour in longitude,
           11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s
           in distance.

         \sstitem
          The original algorithm is expressed in terms of the obsolete
           timescale Ephemeris Time.  Either TDB or TT can be used, but
           not UT without incurring significant errors (30 arcsec at
           the present time) due to the Moon\texttt{'} s 0.5 arcsec/sec movement.

         \sstitem
          The algorithm is based on pre IAU 1976 standards.  However,
           the result has been moved onto the new (FK5) equinox, an
           adjustment which is in any case much smaller than the
           intrinsic accuracy of the procedure.

         \sstitem
          Velocity is obtained by a complete analytical differentiation
           of the Meeus model.
      }
   }
}
\sstroutine{
   palDrange
}{
   Normalize angle into range $+$/- pi
}{
   \sstdescription{
      The result is \texttt{"} angle\texttt{"}  expressed in the range $+$/- pi. If the
      supplied value for \texttt{"} angle\texttt{"}  is equal to $+$/- pi, it is returned
      unchanged.
   }
   \sstinvocation{
      palDrange( double angle )
   }
   \sstarguments{
      \sstsubsection{
         angle = double (Given)
      }{
         The angle in radians.
      }
   }
}
\sstroutine{
   palDt
}{
   Estimate the offset between dynamical time and UT
}{
   \sstdescription{
      Estimate the offset between dynamical time and Universal Time
      for a given historical epoch.
   }
   \sstinvocation{
      double palDt( double epoch );
   }
   \sstarguments{
      \sstsubsection{
         epoch = double (Given)
      }{
         Julian epoch (e.g. 1850.0)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         palDt = double
      }{
         Rough estimate of ET-UT (after 1984, TT-UT) at the
         given epoch, in seconds.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Depending on the epoch, one of three parabolic approximations
           is used:

      }
          before 979    Stephenson \& Morrison\texttt{'} s 390 BC to AD 948 model
          979 to 1708   Stephenson \& Morrison\texttt{'} s 948 to 1600 model
          after 1708    McCarthy \& Babcock\texttt{'} s post-1650 model

        The breakpoints are chosen to ensure continuity:  they occur
        at places where the adjacent models give the same answer as
        each other.
      \sstitemlist{

         \sstitem
          The accuracy is modest, with errors of up to 20 sec during
           the interval since 1650, rising to perhaps 30 min by 1000 BC.
           Comparatively accurate values from AD 1600 are tabulated in
           the Astronomical Almanac (see section K8 of the 1995 AA).

         \sstitem
          The use of double-precision for both argument and result is
           purely for compatibility with other SLALIB time routines.

         \sstitem
          The models used are based on a lunar tidal acceleration value
           of -26.00 arcsec per century.
      }
   }
   \sstdiytopic{
      See Also
   }{
      Explanatory Supplement to the Astronomical Almanac,
      ed P.K.Seidelmann, University Science Books (1992),
      section 2.553, p83.  This contains references to
      the Stephenson \& Morrison and McCarthy \& Babcock
      papers.
   }
}
\sstroutine{
   palDtp2s
}{
   Tangent plane to spherical coordinates
}{
   \sstdescription{
      Transform tangent plane coordinates into spherical.
   }
   \sstinvocation{
      palDtp2s( double xi, double eta, double raz, double decz,
                double $*$ra, double $*$dec);
   }
   \sstarguments{
      \sstsubsection{
         xi = double (Given)
      }{
         First rectangular coordinate on tangent plane (radians)
      }
      \sstsubsection{
         eta = double (Given)
      }{
         Second rectangular coordinate on tangent plane (radians)
      }
      \sstsubsection{
         raz = double (Given)
      }{
         RA spherical coordinate of tangent point (radians)
      }
      \sstsubsection{
         decz = double (Given)
      }{
         Dec spherical coordinate of tangent point (radians)
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         RA spherical coordinate of point to be projected (radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Dec spherical coordinate of point to be projected (radians)
      }
   }
}
\sstroutine{
   palDtps2c
}{
   Determine RA,Dec of tangent point from coordinates
}{
   \sstdescription{
      From the tangent plane coordinates of a star of known RA,Dec,
      determine the RA,Dec of the tangent point.
   }
   \sstinvocation{
      palDtps2c( double xi, double eta, double ra, double dec,
                 double $*$ raz1, double decz1,
                 double $*$ raz2, double decz2, int $*$n);
   }
   \sstarguments{
      \sstsubsection{
         xi = double (Given)
      }{
         First rectangular coordinate on tangent plane (radians)
      }
      \sstsubsection{
         eta = double (Given)
      }{
         Second rectangular coordinate on tangent plane (radians)
      }
      \sstsubsection{
         ra = double (Given)
      }{
         RA spherical coordinate of star (radians)
      }
      \sstsubsection{
         dec = double (Given)
      }{
         Dec spherical coordinate of star (radians)
      }
      \sstsubsection{
         raz1 = double $*$ (Returned)
      }{
         RA spherical coordinate of tangent point, solution 1 (radians)
      }
      \sstsubsection{
         decz1 = double $*$ (Returned)
      }{
         Dec spherical coordinate of tangent point, solution 1 (radians)
      }
      \sstsubsection{
         raz2 = double $*$ (Returned)
      }{
         RA spherical coordinate of tangent point, solution 2 (radians)
      }
      \sstsubsection{
         decz2 = double $*$ (Returned)
      }{
         Dec spherical coordinate of tangent point, solution 2 (radians)
      }
      \sstsubsection{
         n = int $*$ (Returned)
      }{
         number of solutions: 0 = no solutions returned (note 2)
                              1 = only the first solution is useful (note 3)
                              2 = both solutions are useful (note 3)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The RAZ1 and RAZ2 values are returned in the range 0-2pi.

         \sstitem
          Cases where there is no solution can only arise near the poles.
           For example, it is clearly impossible for a star at the pole
           itself to have a non-zero XI value, and hence it is
           meaningless to ask where the tangent point would have to be
           to bring about this combination of XI and DEC.

         \sstitem
          Also near the poles, cases can arise where there are two useful
           solutions.  The argument N indicates whether the second of the
           two solutions returned is useful.  N=1 indicates only one useful
           solution, the usual case;  under these circumstances, the second
           solution corresponds to the \texttt{"} over-the-pole\texttt{"}  case, and this is
           reflected in the values of RAZ2 and DECZ2 which are returned.

         \sstitem
          The DECZ1 and DECZ2 values are returned in the range $+$/-pi, but
           in the usual, non-pole-crossing, case, the range is $+$/-pi/2.

         \sstitem
          This routine is the spherical equivalent of the routine sla\_DTPV2C.
      }
   }
}
\sstroutine{
   palDtt
}{
   Return offset between UTC and TT
}{
   \sstdescription{
      Increment to be applied to Coordinated Universal Time UTC to give
      Terrestrial Time TT (formerly Ephemeris Time ET)
   }
   \sstinvocation{
      dtt = palDtt( double utc );
   }
   \sstarguments{
      \sstsubsection{
         utc = double (Given)
      }{
         UTC date as a modified JD (JD-2400000.5)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         dtt = double
      }{
         TT-UTC in seconds
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Consider a comprehensive upgrade to use the time transformations in SOFA\texttt{'} s time
           cookbook:  http://www.iausofa.org/sofa\_ts\_c.pdf.

         \sstitem
          See eraDat for a description of error conditions when calling this function
           with a time outside of the UTC range. This behaviour differs from slaDtt.
      }
   }
}
\sstroutine{
   palEcleq
}{
   Transform from ecliptic coordinates to J2000.0 equatorial coordinates
}{
   \sstdescription{
      Transform from ecliptic coordinate to J2000.0 equatorial coordinates.
   }
   \sstinvocation{
      void palEcleq ( double dl, double db, double date,
                      double $*$dr, double $*$dd );
   }
   \sstarguments{
      \sstsubsection{
         dl = double (Given)
      }{
         Ecliptic longitude (mean of date, IAU 1980 theory, radians)
      }
      \sstsubsection{
         db = double (Given)
      }{
         Ecliptic latitude (mean of date, IAU 1980 theory, radians)
      }
      \sstsubsection{
         date = double (Given)
      }{
         TT as Modified Julian Date (JD-2400000.5). The difference
         between TT and TDB is of the order of a millisecond or two
         (i.e. about 0.02 arc-seconds).
      }
      \sstsubsection{
         dr = double $*$ (Returned)
      }{
         J2000.0 mean RA (radians)
      }
      \sstsubsection{
         dd = double $*$ (Returned)
      }{
         J2000.0 mean Dec (Radians)
      }
   }
}
\sstroutine{
   palEcmat
}{
   Form the equatorial to ecliptic rotation matrix - IAU 2006
   precession model
}{
   \sstdescription{
      The equatorial to ecliptic rotation matrix is found and returned.
      The matrix is in the sense   V(ecl)  =  RMAT $*$ V(equ);  the
      equator, equinox and ecliptic are mean of date.
   }
   \sstinvocation{
      palEcmat( double date, double rmat[3][3] )
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT as Modified Julian Date (JD-2400000.5). The difference
         between TT and TDB is of the order of a millisecond or two
         (i.e. about 0.02 arc-seconds).
      }
      \sstsubsection{
         rmat = double[3][3] (Returned)
      }{
         Rotation matrix
      }
   }
}
\sstroutine{
   palEl2ue
}{
   Transform conventional elements into \texttt{"} universal\texttt{"}  form
}{
   \sstdescription{
      Transform conventional osculating elements into \texttt{"} universal\texttt{"}  form.
   }
   \sstinvocation{
      void palEl2ue ( double date, int jform, double epoch, double orbinc,
                      double anode, double perih, double aorq, double e,
                      double aorl, double dm, double u[13], int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         Epoch (TT MJD) of osculation (Note 3)
      }
      \sstsubsection{
         jform = int (Given)
      }{
         Element set actually returned (1-3; Note 6)
      }
      \sstsubsection{
         epoch = double (Given)
      }{
         Epoch of elements (TT MJD)
      }
      \sstsubsection{
         orbinc = double (Given)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode = double (Given)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih = double (Given)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq = double (Given)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e = double (Given)
      }{
         eccentricity
      }
      \sstsubsection{
         aorl = double (Given)
      }{
         mean anomaly or longitude (radians, JFORM=1,2 only)
      }
      \sstsubsection{
         dm = double (Given)
      }{
         daily motion (radians, JFORM=1 only)
      }
      \sstsubsection{
         u = double [13] (Returned)
      }{
         Universal orbital elements (Note 1)
         \sstitemlist{

            \sstitem
               (0)  combined mass (M$+$m)

            \sstitem
               (1)  total energy of the orbit (alpha)

            \sstitem
               (2)  reference (osculating) epoch (t0)

            \sstitem
               (3-5)  position at reference epoch (r0)

            \sstitem
               (6-8)  velocity at reference epoch (v0)

            \sstitem
               (9)  heliocentric distance at reference epoch

            \sstitem
               (10)  r0.v0

            \sstitem
               (11)  date (t)

            \sstitem
               (12)  universal eccentric anomaly (psi) of date, approx
         }
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:  0 = OK
         \sstitemlist{

            \sstitem
                   -1 = illegal JFORM

            \sstitem
                   -2 = illegal E

            \sstitem
                   -3 = illegal AORQ

            \sstitem
                   -4 = illegal DM

            \sstitem
                   -5 = numerical error
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The \texttt{"} universal\texttt{"}  elements are those which define the orbit for the
           purposes of the method of universal variables (see reference).
           They consist of the combined mass of the two bodies, an epoch,
           and the position and velocity vectors (arbitrary reference frame)
           at that epoch.  The parameter set used here includes also various
           quantities that can, in fact, be derived from the other
           information.  This approach is taken to avoiding unnecessary
           computation and loss of accuracy.  The supplementary quantities
           are (i) alpha, which is proportional to the total energy of the
           orbit, (ii) the heliocentric distance at epoch, (iii) the
           outwards component of the velocity at the given epoch, (iv) an
           estimate of psi, the \texttt{"} universal eccentric anomaly\texttt{"}  at a given
           date and (v) that date.

         \sstitem
          The companion routine is palUe2pv.  This takes the set of numbers
           that the present routine outputs and uses them to derive the
           object\texttt{'} s position and velocity.  A single prediction requires one
           call to the present routine followed by one call to palUe2pv;
           for convenience, the two calls are packaged as the routine
           palPlanel.  Multiple predictions may be made by again calling the
           present routine once, but then calling palUe2pv multiple times,
           which is faster than multiple calls to palPlanel.

         \sstitem
          DATE is the epoch of osculation.  It is in the TT timescale
           (formerly Ephemeris Time, ET) and is a Modified Julian Date
           (JD-2400000.5).

         \sstitem
          The supplied orbital elements are with respect to the J2000
           ecliptic and equinox.  The position and velocity parameters
           returned in the array U are with respect to the mean equator and
           equinox of epoch J2000, and are for the perihelion prior to the
           specified epoch.

         \sstitem
          The universal elements returned in the array U are in canonical
           units (solar masses, AU and canonical days).

         \sstitem
          Three different element-format options are available:

      }
        Option JFORM=1, suitable for the major planets:

        EPOCH  = epoch of elements (TT MJD)
        ORBINC = inclination i (radians)
        ANODE  = longitude of the ascending node, big omega (radians)
        PERIH  = longitude of perihelion, curly pi (radians)
        AORQ   = mean distance, a (AU)
        E      = eccentricity, e (range 0 to $<$1)
        AORL   = mean longitude L (radians)
        DM     = daily motion (radians)

        Option JFORM=2, suitable for minor planets:

        EPOCH  = epoch of elements (TT MJD)
        ORBINC = inclination i (radians)
        ANODE  = longitude of the ascending node, big omega (radians)
        PERIH  = argument of perihelion, little omega (radians)
        AORQ   = mean distance, a (AU)
        E      = eccentricity, e (range 0 to $<$1)
        AORL   = mean anomaly M (radians)

        Option JFORM=3, suitable for comets:

        EPOCH  = epoch of perihelion (TT MJD)
        ORBINC = inclination i (radians)
        ANODE  = longitude of the ascending node, big omega (radians)
        PERIH  = argument of perihelion, little omega (radians)
        AORQ   = perihelion distance, q (AU)
        E      = eccentricity, e (range 0 to 10)

      \sstitemlist{

         \sstitem
          Unused elements (DM for JFORM=2, AORL and DM for JFORM=3) are
           not accessed.

         \sstitem
          The algorithm was originally adapted from the EPHSLA program of
           D.H.P.Jones (private communication, 1996).  The method is based
           on Stumpff\texttt{'} s Universal Variables.
      }
   }
   \sstdiytopic{
      See Also
   }{
      Everhart \& Pitkin, Am.J.Phys. 51, 712 (1983).
   }
}
\sstroutine{
   palEpco
}{
   Convert an epoch into the appropriate form - \texttt{'} B\texttt{'}  or \texttt{'} J\texttt{'}
}{
   \sstdescription{
      Converts a Besselian or Julian epoch to a Julian or Besselian
      epoch.
   }
   \sstinvocation{
      double palEpco( char k0, char k, double e );
   }
   \sstarguments{
      \sstsubsection{
         k0 = char (Given)
      }{
         Form of result: \texttt{'} B\texttt{'} =Besselian, \texttt{'} J\texttt{'} =Julian
      }
      \sstsubsection{
         k = char (Given)
      }{
         Form of given epoch: \texttt{'} B\texttt{'}  or \texttt{'} J\texttt{'} .
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The result is always either equal to or very close to
           the given epoch E.  The routine is required only in
           applications where punctilious treatment of heterogeneous
           mixtures of star positions is necessary.

         \sstitem
          k and k0 are case insensitive. This differes slightly from the
           Fortran SLA implementation.

         \sstitem
          k and k0 are not validated. They are interpreted as follows:
           o If k0 and k are the same the result is e
           o If k0 is \texttt{'} b\texttt{'}  or \texttt{'} B\texttt{'}  and k isn\texttt{'} t the conversion is J to B.
           o In all other cases, the conversion is B to J.
      }
   }
}
\sstroutine{
   palEpv
}{
   Earth position and velocity with respect to the BCRS
}{
   \sstdescription{
      Earth position and velocity, heliocentric and barycentric, with
      respect to the Barycentric Celestial Reference System.
   }
   \sstinvocation{
      void palEpv( double date, double ph[3], double vh[3],
                   double pb[3], double vb[3] );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         Date, TDB Modified Julian Date (JD-2400000.5)
      }
      \sstsubsection{
         ph = double [3] (Returned)
      }{
         Heliocentric Earth position (AU)
      }
      \sstsubsection{
         vh = double [3] (Returned)
      }{
         Heliocentric Earth velocity (AU/day)
      }
      \sstsubsection{
         pb = double [3] (Returned)
      }{
         Barycentric Earth position (AU)
      }
      \sstsubsection{
         vb = double [3] (Returned)
      }{
         Barycentric Earth velocity (AU/day)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          See eraEpv00 for details on accuracy

         \sstitem
          Note that the status argument from eraEpv00 is ignored
      }
   }
}
\sstroutine{
   palEtrms
}{
   Compute the E-terms vector
}{
   \sstdescription{
      Computes the E-terms (elliptic component of annual aberration)
      vector.

      Note the use of the J2000 aberration constant (20.49552 arcsec).
      This is a reflection of the fact that the E-terms embodied in
      existing star catalogues were computed from a variety of
      aberration constants.  Rather than adopting one of the old
      constants the latest value is used here.
   }
   \sstinvocation{
      void palEtrms ( double ep, double ev[3] );
   }
   \sstarguments{
      \sstsubsection{
         ep = double (Given)
      }{
         Besselian epoch
      }
      \sstsubsection{
         ev = double [3] (Returned)
      }{
         E-terms as (dx,dy,dz)
      }
   }
   \sstdiytopic{
      See also
   }{
      \sstitemlist{

         \sstitem
          Smith, C.A. et al., 1989.  Astr.J. 97, 265.

         \sstitem
          Yallop, B.D. et al., 1989.  Astr.J. 97, 274.
      }
   }
}
\sstroutine{
   palEqecl
}{
   Transform from J2000.0 equatorial coordinates to ecliptic coordinates
}{
   \sstdescription{
      Transform from J2000.0 equatorial coordinates to ecliptic coordinates.
   }
   \sstinvocation{
      void palEqecl( double dr, double dd, double date,
                     double $*$dl, double $*$db);
   }
   \sstarguments{
      \sstsubsection{
         dr = double (Given)
      }{
         J2000.0 mean RA (radians)
      }
      \sstsubsection{
         dd = double (Given)
      }{
         J2000.0 mean Dec (Radians)
      }
      \sstsubsection{
         date = double (Given)
      }{
         TT as Modified Julian Date (JD-2400000.5). The difference
         between TT and TDB is of the order of a millisecond or two
         (i.e. about 0.02 arc-seconds).
      }
      \sstsubsection{
         dl = double $*$ (Returned)
      }{
         Ecliptic longitude (mean of date, IAU 1980 theory, radians)
      }
      \sstsubsection{
         db = double $*$ (Returned)
      }{
         Ecliptic latitude (mean of date, IAU 1980 theory, radians)
      }
   }
}
\sstroutine{
   palEqgal
}{
   Convert from J2000.0 equatorial coordinates to Galactic
}{
   \sstdescription{
      Transformation from J2000.0 equatorial coordinates
      to IAU 1958 galactic coordinates.
   }
   \sstinvocation{
      void palEqgal ( double dr, double dd, double $*$dl, double $*$db );
   }
   \sstarguments{
      \sstsubsection{
         dr = double (Given)
      }{
         J2000.0 RA (radians)
      }
      \sstsubsection{
         dd = double (Given)
      }{
         J2000.0 Dec (radians
      }
      \sstsubsection{
         dl = double $*$ (Returned)
      }{
         Galactic longitude (radians).
      }
      \sstsubsection{
         db = double $*$ (Returned)
      }{
         Galactic latitude (radians).
      }
   }
   \sstnotes{
      The equatorial coordinates are J2000.0.  Use the routine
      palGe50 if conversion to B1950.0 \texttt{'} FK4\texttt{'}  coordinates is
      required.
   }
   \sstdiytopic{
      See Also
   }{
      Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960)
   }
}
\sstroutine{
   palEvp
}{
   Returns the barycentric and heliocentric velocity and position of the
   Earth
}{
   \sstdescription{
      Returns the barycentric and heliocentric velocity and position of the
      Earth at a given epoch, given with respect to a specified equinox.
      For information about accuracy, see the function eraEpv00.
   }
   \sstinvocation{
      void palEvp( double date, double deqx, double dvb[3], double dpb[3],
                   double dvh[3], double dph[3] )
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TDB (loosely ET) as a Modified Julian Date (JD-2400000.5)
      }
      \sstsubsection{
         deqx = double (Given)
      }{
         Julian epoch (e.g. 2000.0) of mean equator and equinox of the
         vectors returned.  If deqx $<$= 0.0, all vectors are referred to the
         mean equator and equinox (FK5) of epoch date.
      }
      \sstsubsection{
         dvb = double[3] (Returned)
      }{
         Barycentric velocity (AU/s, AU)
      }
      \sstsubsection{
         dpb = double[3] (Returned)
      }{
         Barycentric position (AU/s, AU)
      }
      \sstsubsection{
         dvh = double[3] (Returned)
      }{
         heliocentric velocity (AU/s, AU)
      }
      \sstsubsection{
         dph = double[3] (Returned)
      }{
         Heliocentric position (AU/s, AU)
      }
   }
}
\sstroutine{
   palFitxy
}{
   Fit a linear model to relate two sets of [x,y] coordinates
}{
   \sstdescription{
      Fits a linear model to relate two sets of [X,Y] coordinates.
   }
   \sstinvocation{
      palFitxy ( int itype, int np, double xye[][2], double xym[][2],
                 double coeffs[6], int $*$j )
   }
   \sstarguments{
      \sstsubsection{
         itype = int (Given)
      }{
         type of model: 4 or 6 (note 1)
      }
      \sstsubsection{
         np = int (Given)
      }{
         number of samples (note 2)
      }
      \sstsubsection{
         xye = double[np][2] (Given)
      }{
         expected [x,y] for each sample
      }
      \sstsubsection{
         xym = double[np][2] (Given)
      }{
         measured [x,y] for each sample
      }
      \sstsubsection{
         coeffs = double[6] (Returned)
      }{
         coefficients of model (note 3)
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         status:
         \sstitemlist{

            \sstitem
             0 = OK

            \sstitem
             -1 = illegal itype

            \sstitem
             -2 = insufficient data

            \sstitem
             -3 = no solution
         }
      }
   }
   \sstnotes{
      1)  itype, which must be either 4 or 6, selects the type of model
          fitted.  Both allowed itype values produce a model coeffs which
          consists of six coefficients, namely the zero points and, for
          each of xe and ye, the coefficient of xm and ym.  For itype=6,
          all six coefficients are independent, modelling squash and shear
          as well as origin, scale, and orientation.  However, itype=4
          selects the \texttt{"} solid body rotation\texttt{"}  option;  the model coeffs
          still consists of the same six coefficients, but now two of
          them are used twice (appropriately signed).  Origin, scale
          and orientation are still modelled, but not squash or shear -
          the units of x and y have to be the same.

      2)  For itype=4, np must be at least 2.  For itype=6, np must be at
          least 3.

      3)  The model is returned in the array coeffs.  Naming the
          elements of coeffs as follows:
      \newline
      \newline
      \hspace*{1.5 em}
             coeffs[0] = A
      \newline
      \hspace*{1.5 em}
             coeffs[1] = B
      \newline
      \hspace*{1.5 em}
             coeffs[2] = C
      \newline
      \hspace*{1.5 em}
             coeffs[3] = D
      \newline
      \hspace*{1.5 em}
             coeffs[4] = E
      \newline
      \hspace*{1.5 em}
             coeffs[5] = F
      \newline
      \newline
          the model is:
      \newline
      \newline
      \hspace*{1.5 em}
             xe = A $+$ B $*$ xm $+$ C $*$ ym
      \newline
      \hspace*{1.5 em}
             ye = D $+$ E $*$ xm $+$ F $*$ ym
      \newline
      \newline
          For the \texttt{"} solid body rotation\texttt{"}  option (itype=4), the
          magnitudes of B and F, and of C and E, are equal.  The
          signs of these coefficients depend on whether there is a
          sign reversal between xe,ye and xm,ym;  fits are performed
          with and without a sign reversal and the best one chosen.

      4)  Error status values j=-1 and -2 leave coeffs unchanged;
          if j=-3 coeffs may have been changed.
   }
   \sstdiytopic{
      See also
   }{
      palPxy, palInvf, palXy2xy and palDcmpf
   }
}
\sstroutine{
   palFk45z
}{
   Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
   proper motion in the FK5 frame
}{
   \sstdescription{
      Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero
      proper motion in the FK5 frame (double precision)

      This function converts stars from the Bessel-Newcomb, FK4
      system to the IAU 1976, FK5, Fricke system, in such a
      way that the FK5 proper motion is zero.  Because such a star
      has, in general, a non-zero proper motion in the FK4 system,
      the routine requires the epoch at which the position in the
      FK4 system was determined.

      The method is from Appendix 2 of Ref 1, but using the constants
      of Ref 4.
   }
   \sstinvocation{
      palFk45z( double r1950, double d1950, double bepoch, double $*$r2000,
                double $*$d2000 )
   }
   \sstarguments{
      \sstsubsection{
         r1950 = double (Given)
      }{
         B1950.0 FK4 RA at epoch (radians).
      }
      \sstsubsection{
         d1950 = double (Given)
      }{
         B1950.0 FK4 Dec at epoch (radians).
      }
      \sstsubsection{
         bepoch = double (Given)
      }{
         Besselian epoch (e.g. 1979.3)
      }
      \sstsubsection{
         r2000 = double (Returned)
      }{
         J2000.0 FK5 RA (Radians).
      }
      \sstsubsection{
         d2000 = double (Returned)
      }{
         J2000.0 FK5 Dec(Radians).
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The epoch BEPOCH is strictly speaking Besselian, but if a
         Julian epoch is supplied the result will be affected only to
         a negligible extent.

         \sstitem
          Conversion from Besselian epoch 1950.0 to Julian epoch 2000.0
         only is provided for.  Conversions involving other epochs will
         require use of the appropriate precession, proper motion, and
         E-terms routines before and/or after palFk45z is called.

         \sstitem
          In the FK4 catalogue the proper motions of stars within 10
         degrees of the poles do not embody the differential E-term effect
         and should, strictly speaking, be handled in a different manner
         from stars outside these regions. However, given the general lack
         of homogeneity of the star data available for routine astrometry,
         the difficulties of handling positions that may have been
         determined from astrometric fields spanning the polar and non-polar
         regions, the likelihood that the differential E-terms effect was not
         taken into account when allowing for proper motion in past
         astrometry, and the undesirability of a discontinuity in the
         algorithm, the decision has been made in this routine to include the
         effect of differential E-terms on the proper motions for all stars,
         whether polar or not.  At epoch 2000, and measuring on the sky rather
         than in terms of dRA, the errors resulting from this simplification
         are less than 1 milliarcsecond in position and 1 milliarcsecond per
         century in proper motion.
      }
   }
   \sstdiytopic{
      References
   }{
      \sstitemlist{

         \sstitem
          Aoki,S., et al, 1983.  Astron.Astrophys., 128, 263.

         \sstitem
          Smith, C.A. et al, 1989.  \texttt{"} The transformation of astrometric
           catalog systems to the equinox J2000.0\texttt{"} .  Astron.J. 97, 265.

         \sstitem
          Yallop, B.D. et al, 1989.  \texttt{"} Transformation of mean star places
           from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space\texttt{"} .
           Astron.J. 97, 274.

         \sstitem
          Seidelmann, P.K. (ed), 1992.  \texttt{"} Explanatory Supplement to
           the Astronomical Almanac\texttt{"} , ISBN 0-935702-68-7.
      }
   }
}
\sstroutine{
   palFk524
}{
   Convert J2000.0 FK5 star data to B1950.0 FK4
}{
   \sstdescription{
      This function converts stars from the IAU 1976, FK5, Fricke
      system, to the Bessel-Newcomb, FK4 system.  The precepts
      of Smith et al (Ref 1) are followed, using the implementation
      by Yallop et al (Ref 2) of a matrix method due to Standish.
      Kinoshita\texttt{'} s development of Andoyer\texttt{'} s post-Newcomb precession is
      used.  The numerical constants from Seidelmann et al (Ref 3) are
      used canonically.
   }
   \sstinvocation{
      palFk524( double r2000, double d2000, double dr2000, double dd2000,
                double p2000, double v2000, double $*$r1950, double $*$d1950,
                double $*$dr1950, double $*$dd1950, double $*$p1950, double $*$v1950 )
   }
   \sstarguments{
      \sstsubsection{
         r2000 = double (Given)
      }{
         J2000.0 FK5 RA (radians).
      }
      \sstsubsection{
         d2000 = double (Given)
      }{
         J2000.0 FK5 Dec (radians).
      }
      \sstsubsection{
         dr2000 = double (Given)
      }{
         J2000.0 FK5 RA proper motion (rad/Jul.yr)
      }
      \sstsubsection{
         dd2000 = double (Given)
      }{
         J2000.0 FK5 Dec proper motion (rad/Jul.yr)
      }
      \sstsubsection{
         p2000 = double (Given)
      }{
         J2000.0 FK5 parallax (arcsec)
      }
      \sstsubsection{
         v2000 = double (Given)
      }{
         J2000.0 FK5 radial velocity (km/s, $+$ve = moving away)
      }
      \sstsubsection{
         r1950 = double $*$ (Returned)
      }{
         B1950.0 FK4 RA (radians).
      }
      \sstsubsection{
         d1950 = double $*$ (Returned)
      }{
         B1950.0 FK4 Dec (radians).
      }
      \sstsubsection{
         dr1950 = double $*$ (Returned)
      }{
         B1950.0 FK4 RA proper motion (rad/Jul.yr)
      }
      \sstsubsection{
         dd1950 = double $*$ (Returned)
      }{
         B1950.0 FK4 Dec proper motion (rad/Jul.yr)
      }
      \sstsubsection{
         p1950 = double $*$ (Returned)
      }{
         B1950.0 FK4 parallax (arcsec)
      }
      \sstsubsection{
         v1950 = double $*$ (Returned)
      }{
         B1950.0 FK4 radial velocity (km/s, $+$ve = moving away)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The proper motions in RA are dRA/dt rather than
         cos(Dec)$*$dRA/dt, and are per year rather than per century.

         \sstitem
          Note that conversion from Julian epoch 2000.0 to Besselian
         epoch 1950.0 only is provided for.  Conversions involving
         other epochs will require use of the appropriate precession,
         proper motion, and E-terms routines before and/or after
         FK524 is called.

         \sstitem
          In the FK4 catalogue the proper motions of stars within
         10 degrees of the poles do not embody the differential
         E-term effect and should, strictly speaking, be handled
         in a different manner from stars outside these regions.
         However, given the general lack of homogeneity of the star
         data available for routine astrometry, the difficulties of
         handling positions that may have been determined from
         astrometric fields spanning the polar and non-polar regions,
         the likelihood that the differential E-terms effect was not
         taken into account when allowing for proper motion in past
         astrometry, and the undesirability of a discontinuity in
         the algorithm, the decision has been made in this routine to
         include the effect of differential E-terms on the proper
         motions for all stars, whether polar or not.  At epoch 2000,
         and measuring on the sky rather than in terms of dRA, the
         errors resulting from this simplification are less than
         1 milliarcsecond in position and 1 milliarcsecond per
         century in proper motion.
      }
   }
   \sstdiytopic{
      References
   }{
      \sstitemlist{

         \sstitem
          Smith, C.A. et al, 1989.  \texttt{"} The transformation of astrometric
           catalog systems to the equinox J2000.0\texttt{"} .  Astron.J. 97, 265.

         \sstitem
          Yallop, B.D. et al, 1989.  \texttt{"} Transformation of mean star places
           from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space\texttt{"} .
           Astron.J. 97, 274.

         \sstitem
          Seidelmann, P.K. (ed), 1992.  \texttt{"} Explanatory Supplement to
           the Astronomical Almanac\texttt{"} , ISBN 0-935702-68-7.
      }
   }
}
\sstroutine{
   palFk54z
}{
   Convert a J2000.0 FK5 star position to B1950.0 FK4 assuming
   zero proper motion and parallax
}{
   \sstdescription{
      This function converts star positions from the IAU 1976,
      FK5, Fricke system to the Bessel-Newcomb, FK4 system.
   }
   \sstinvocation{
      palFk54z( double r2000, double d2000, double bepoch, double $*$r1950,
                double $*$d1950, double $*$dr1950, double $*$dd1950 )
   }
   \sstarguments{
      \sstsubsection{
         r2000 = double (Given)
      }{
         J2000.0 FK5 RA (radians).
      }
      \sstsubsection{
         d2000 = double (Given)
      }{
         J2000.0 FK5 Dec (radians).
      }
      \sstsubsection{
         bepoch = double (Given)
      }{
         Besselian epoch (e.g. 1950.0).
      }
      \sstsubsection{
         r1950 = double $*$ (Returned)
      }{
         B1950 FK4 RA (radians) at epoch \texttt{"} bepoch\texttt{"} .
      }
      \sstsubsection{
         d1950 = double $*$ (Returned)
      }{
         B1950 FK4 Dec (radians) at epoch \texttt{"} bepoch\texttt{"} .
      }
      \sstsubsection{
         dr1950 = double $*$ (Returned)
      }{
         B1950 FK4 proper motion (RA) (radians/trop.yr)).
      }
      \sstsubsection{
         dr1950 = double $*$ (Returned)
      }{
         B1950 FK4 proper motion (Dec) (radians/trop.yr)).
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The proper motion in RA is dRA/dt rather than cos(Dec)$*$dRA/dt.

         \sstitem
          Conversion from Julian epoch 2000.0 to Besselian epoch 1950.0
         only is provided for.  Conversions involving other epochs will
         require use of the appropriate precession functions before and
         after this function is called.

         \sstitem
          The FK5 proper motions, the parallax and the radial velocity
          are presumed zero.

         \sstitem
          It is the intention that FK5 should be a close approximation
         to an inertial frame, so that distant objects have zero proper
         motion;  such objects have (in general) non-zero proper motion
         in FK4, and this function returns those fictitious proper
         motions.

         \sstitem
          The position returned by this function is in the B1950
         reference frame but at Besselian epoch BEPOCH.  For comparison
         with catalogues the \texttt{"} bepoch\texttt{"}  argument will frequently be 1950.0.
      }
   }
}
\sstroutine{
   palGaleq
}{
   Convert from galactic to J2000.0 equatorial coordinates
}{
   \sstdescription{
      Transformation from IAU 1958 galactic coordinates to
      J2000.0 equatorial coordinates.
   }
   \sstinvocation{
      void palGaleq ( double dl, double db, double $*$dr, double $*$dd );
   }
   \sstarguments{
      \sstsubsection{
         dl = double (Given)
      }{
         Galactic longitude (radians).
      }
      \sstsubsection{
         db = double (Given)
      }{
         Galactic latitude (radians).
      }
      \sstsubsection{
         dr = double $*$ (Returned)
      }{
         J2000.0 RA (radians)
      }
      \sstsubsection{
         dd = double $*$ (Returned)
      }{
         J2000.0 Dec (radians)
      }
   }
   \sstnotes{
      The equatorial coordinates are J2000.0.  Use the routine
      palGe50 if conversion to B1950.0 \texttt{'} FK4\texttt{'}  coordinates is
      required.
   }
   \sstdiytopic{
      See Also
   }{
      Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960)
   }
}
\sstroutine{
   palGalsup
}{
   Convert from galactic to supergalactic coordinates
}{
   \sstdescription{
      Transformation from IAU 1958 galactic coordinates to
      de Vaucouleurs supergalactic coordinates.
   }
   \sstinvocation{
      void palGalsup ( double dl, double db, double $*$dsl, double $*$dsb );
   }
   \sstarguments{
      \sstsubsection{
         dl = double (Given)
      }{
         Galactic longitude.
      }
      \sstsubsection{
         db = double (Given)
      }{
         Galactic latitude.
      }
      \sstsubsection{
         dsl = double $*$ (Returned)
      }{
         Supergalactic longitude.
      }
      \sstsubsection{
         dsb = double $*$ (Returned)
      }{
         Supergalactic latitude.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          de Vaucouleurs, de Vaucouleurs, \& Corwin, Second Reference
           Catalogue of Bright Galaxies, U. Texas, page 8.

         \sstitem
          Systems \& Applied Sciences Corp., Documentation for the
           machine-readable version of the above catalogue,
           Contract NAS 5-26490.

      }
      (These two references give different values for the galactic
      longitude of the supergalactic origin.  Both are wrong;  the
      correct value is L2=137.37.)
   }
}
\sstroutine{
   palGe50
}{
   Transform Galactic Coordinate to B1950 FK4
}{
   \sstdescription{
      Transformation from IAU 1958 galactic coordinates to
      B1950.0 \texttt{'} FK4\texttt{'}  equatorial coordinates.
   }
   \sstinvocation{
      palGe50( double dl, double db, double $*$dr, double $*$dd );
   }
   \sstarguments{
      \sstsubsection{
         dl = double (Given)
      }{
         Galactic longitude (radians)
      }
      \sstsubsection{
         db = double (Given)
      }{
         Galactic latitude (radians)
      }
      \sstsubsection{
         dr = double $*$ (Returned)
      }{
         B9150.0 FK4 RA.
      }
      \sstsubsection{
         dd = double $*$ (Returned)
      }{
         B1950.0 FK4 Dec.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The equatorial coordinates are B1950.0 \texttt{'} FK4\texttt{'} . Use the routine
         palGaleq if conversion to J2000.0 coordinates is required.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960)
      }
   }
}
\sstroutine{
   palGeoc
}{
   Convert geodetic position to geocentric
}{
   \sstdescription{
      Convert geodetic position to geocentric.
   }
   \sstinvocation{
      void palGeoc( double p, double h, double $*$ r, double $*$z );
   }
   \sstarguments{
      \sstsubsection{
         p = double (Given)
      }{
         latitude (radians)
      }
      \sstsubsection{
         h = double (Given)
      }{
         height above reference spheroid (geodetic, metres)
      }
      \sstsubsection{
         r = double $*$ (Returned)
      }{
         distance from Earth axis (AU)
      }
      \sstsubsection{
         z = double $*$ (Returned)
      }{
         distance from plane of Earth equator (AU)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Geocentric latitude can be obtained by evaluating atan2(z,r)

         \sstitem
          Uses WGS84 reference ellipsoid and calls eraGd2gc
      }
   }
}
\sstroutine{
   palIntin
}{
   Convert free-format input into an integer
}{
   \sstdescription{
      Extracts a number from an input string starting at the specified
      index.
   }
   \sstinvocation{
      void palIntin( const char $*$ string, int $*$nstrt,
                      long $*$ireslt, int $*$jflag );
   }
   \sstarguments{
      \sstsubsection{
         string = const char $*$ (Given)
      }{
         String containing number to be decoded.
      }
      \sstsubsection{
         nstrt = int $*$ (Given and Returned)
      }{
         Character number indicating where decoding should start.
         On output its value is updated to be the location of the
         possible next value. For compatibility with SLA the first
         character is index 1.
      }
      \sstsubsection{
         ireslt = long $*$ (Returned)
      }{
         Result. Not updated when jflag=1.
      }
      \sstsubsection{
         jflag = int $*$ (Returned)
      }{
         status: -1 = -OK, 0 = $+$OK, 1 = null, 2 = error
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses the strtol() system call to do the parsing. This may lead to
           subtle differences when compared to the SLA/F parsing.

         \sstitem
          Commas are recognized as a special case and are skipped if one happens
           to be the next character when updating nstrt. Additionally the output
           nstrt position will skip past any trailing space.

         \sstitem
          If no number can be found flag will be set to 1.

         \sstitem
          If the number overflows or underflows jflag will be set to 2. For overflow
           the returned result will have the value LONG\_MAX, for underflow it
           will have the value LONG\_MIN.
      }
   }
}
\sstroutine{
   palInvf
}{
   Invert a linear model of the type produced by the palFitxy routine
}{
   \sstdescription{
      The models relate two sets of [x,y] coordinates as follows.
      Naming the elements of fwds:
      \newline
      \newline
      \hspace*{1.5 em}
         fwds[0] = A
      \newline
      \hspace*{1.5 em}
         fwds[1] = B
      \newline
      \hspace*{1.5 em}
         fwds[2] = C
      \newline
      \hspace*{1.5 em}
         fwds[3] = D
      \newline
      \hspace*{1.5 em}
         fwds[4] = E
      \newline
      \hspace*{1.5 em}
         fwds[5] = F
      \newline
      \newline
      where two sets of coordinates [x1,y1] and [x2,y2] are related
      thus:
      \newline
      \newline
      \hspace*{1.5 em}
         x2 = A $+$ B $*$ x1 $+$ C $*$ y1
      \newline
      \hspace*{1.5 em}
         y2 = D $+$ E $*$ x1 $+$ F $*$ y1
      \newline
      \newline
      the present routine generates a new set of coefficients:
      \newline
      \newline
      \hspace*{1.5 em}
         bkwds[0] = P
      \newline
      \hspace*{1.5 em}
         bkwds[1] = Q
      \newline
      \hspace*{1.5 em}
         bkwds[2] = R
      \newline
      \hspace*{1.5 em}
         bkwds[3] = S
      \newline
      \hspace*{1.5 em}
         bkwds[4] = T
      \newline
      \hspace*{1.5 em}
         bkwds[5] = U
      \newline
      \newline
      such that:
      \newline
      \newline
      \hspace*{1.5 em}
         x1 = P $+$ Q $*$ x2 $+$ R $*$ y2
      \newline
      \hspace*{1.5 em}
         y1 = S $+$ T $*$ x2 $+$ U $*$ y2
      \newline
      \newline
      Two successive calls to palInvf will thus deliver a set
      of coefficients equal to the starting values.
   }
   \sstinvocation{
      palInvf ( double fwds[6], double bkwds[6], int $*$j )
   }
   \sstarguments{
      \sstsubsection{
         fwds = double[6] (Given)
      }{
         model coefficients
      }
      \sstsubsection{
         bkwds = double[6] (Returned)
      }{
         inverse model
      }
      \sstsubsection{
         j = int (Returned)
      }{
         status: 0 = OK, -1 = no inverse
      }
   }
   \sstdiytopic{
      See also
   }{
      palFitxy, palPxy, palXy2xy and palDcmpf
   }
}
\sstroutine{
   palMap
}{
   Convert star RA,Dec from mean place to geocentric apparent
}{
   \sstdescription{
      Convert star RA,Dec from mean place to geocentric apparent.
   }
   \sstinvocation{
      void palMap( double rm, double dm, double pr, double pd,
                   double px, double rv, double eq, double date,
                   double $*$ra, double $*$da );
   }
   \sstarguments{
      \sstsubsection{
         rm = double (Given)
      }{
         Mean RA (radians)
      }
      \sstsubsection{
         dm = double (Given)
      }{
         Mean declination (radians)
      }
      \sstsubsection{
         pr = double (Given)
      }{
         RA proper motion, changes per Julian year (radians)
      }
      \sstsubsection{
         pd = double (Given)
      }{
         Dec proper motion, changes per Julian year (radians)
      }
      \sstsubsection{
         px = double (Given)
      }{
         Parallax (arcsec)
      }
      \sstsubsection{
         rv = double (Given)
      }{
         Radial velocity (km/s, $+$ve if receding)
      }
      \sstsubsection{
         eq = double (Given)
      }{
         Epoch and equinox of star data (Julian)
      }
      \sstsubsection{
         date = double (Given)
      }{
         TDB for apparent place (JD-2400000.5)
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         Apparent RA (radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Apparent dec (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Calls palMappa and palMapqk

         \sstitem
          The reference systems and timescales used are IAU 2006.

         \sstitem
          EQ is the Julian epoch specifying both the reference frame and
           the epoch of the position - usually 2000.  For positions where
           the epoch and equinox are different, use the routine palPm to
           apply proper motion corrections before using this routine.

         \sstitem
          The distinction between the required TDB and TT is always
           negligible.  Moreover, for all but the most critical
           applications UTC is adequate.

         \sstitem
          The proper motions in RA are dRA/dt rather than cos(Dec)$*$dRA/dt.

         \sstitem
          This routine may be wasteful for some applications because it
           recomputes the Earth position/velocity and the precession-
           nutation matrix each time, and because it allows for parallax
           and proper motion.  Where multiple transformations are to be
           carried out for one epoch, a faster method is to call the
           palMappa routine once and then either the palMapqk routine
           (which includes parallax and proper motion) or palMapqkz (which
           assumes zero parallax and proper motion).

         \sstitem
          The accuracy is sub-milliarcsecond, limited by the
           precession-nutation model (see palPrenut for details).

         \sstitem
          The accuracy is further limited by the routine palEvp, called
           by palMappa, which computes the Earth position and velocity.
           See eraEpv00 for details on that calculation.
      }
   }
}
\sstroutine{
   palMappa
}{
   Compute parameters needed by palAmpqk and palMapqk
}{
   \sstdescription{
      Compute star-independent parameters in preparation for
      transformations between mean place and geocentric apparent place.

      The parameters produced by this function are required in the
      parallax, aberration, and nutation/bias/precession parts of the
      mean/apparent transformations.

      The reference systems and timescales used are IAU 2006.
   }
   \sstinvocation{
      void palMappa( double eq, double date, double amprms[21] )
   }
   \sstarguments{
      \sstsubsection{
         eq = double (Given)
      }{
         epoch of mean equinox to be used (Julian)
      }
      \sstsubsection{
         date = double (Given)
      }{
         TDB (JD-2400000.5)
      }
      \sstsubsection{
         amprms =   double[21]  (Returned)
      }{
         star-independent mean-to-apparent parameters:
         \sstitemlist{

            \sstitem
             (0)      time interval for proper motion (Julian years)

            \sstitem
             (1-3)    barycentric position of the Earth (AU)

            \sstitem
             (4-6)    heliocentric direction of the Earth (unit vector)

            \sstitem
             (7)      (Schwarzschild radius of Sun)/(Sun-Earth distance)

            \sstitem
             (8-10)   abv: barycentric Earth velocity in units of c

            \sstitem
             (11)     sqrt(1-v$\wedge$2) where v=modulus(abv)

            \sstitem
             (12-20)  precession/nutation (3,3) matrix
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          For date, the distinction between the required TDB and TT
         is always negligible.  Moreover, for all but the most
         critical applications UTC is adequate.

         \sstitem
          The vector amprms(1-3) is referred to the mean equinox and
         equator of epoch eq.

         \sstitem
          The parameters amprms produced by this function are used by
         palAmpqk, palMapqk and palMapqkz.
      }
   }
}
\sstroutine{
   palMapqk
}{
   Quick mean to apparent place
}{
   \sstdescription{
      Quick mean to apparent place:  transform a star RA,Dec from
      mean place to geocentric apparent place, given the
      star-independent parameters.

      Use of this routine is appropriate when efficiency is important
      and where many star positions, all referred to the same equator
      and equinox, are to be transformed for one epoch.  The
      star-independent parameters can be obtained by calling the
      palMappa routine.

      If the parallax and proper motions are zero the palMapqkz
      routine can be used instead.
   }
   \sstinvocation{
      void palMapqk ( double rm, double dm, double pr, double pd,
                      double px, double rv, double amprms[21],
                      double $*$ra, double $*$da );
   }
   \sstarguments{
      \sstsubsection{
         rm = double (Given)
      }{
         Mean RA (radians)
      }
      \sstsubsection{
         dm = double (Given)
      }{
         Mean declination (radians)
      }
      \sstsubsection{
         pr = double (Given)
      }{
         RA proper motion, changes per Julian year (radians)
      }
      \sstsubsection{
         pd = double (Given)
      }{
         Dec proper motion, changes per Julian year (radians)
      }
      \sstsubsection{
         px = double (Given)
      }{
         Parallax (arcsec)
      }
      \sstsubsection{
         rv = double (Given)
      }{
         Radial velocity (km/s, $+$ve if receding)
      }
      \sstsubsection{
         amprms = double [21] (Given)
      }{
         Star-independent mean-to-apparent parameters (see palMappa).
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         Apparent RA (radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Apparent dec (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The reference frames and timescales used are post IAU 2006.

         \sstitem
          The mean place rm, dm and the vectors amprms[1-3] and amprms[4-6]
           are referred to the mean equinox and equator of the epoch
           specified when generating the precession/nutation matrix
           amprms[12-20].  In the call to palMappa (q.v.) normally used
           to populate amprms, this epoch is the first argument (eq).

         \sstitem
          Strictly speaking, the routine is not valid for solar-system
           sources, though the error will usually be extremely small.
           However, to prevent gross errors in the case where the
           position of the Sun is specified, the gravitational
           deflection term is restrained within about 920 arcsec of the
           centre of the Sun\texttt{'} s disc.  The term has a maximum value of
           about 1.85 arcsec at this radius, and decreases to zero as
           the centre of the disc is approached.
      }
   }
}
\sstroutine{
   palMapqkz
}{
   Quick mean to apparent place (no proper motion or parallax)
}{
   \sstdescription{
      Quick mean to apparent place:  transform a star RA,dec from
      mean place to geocentric apparent place, given the
      star-independent parameters, and assuming zero parallax
      and proper motion.

      Use of this function is appropriate when efficiency is important
      and where many star positions, all with parallax and proper
      motion either zero or already allowed for, and all referred to
      the same equator and equinox, are to be transformed for one
      epoch.  The star-independent parameters can be obtained by
      calling the palMappa function.

      The corresponding function for the case of non-zero parallax
      and proper motion is palMapqk.
   }
   \sstinvocation{
      void palMapqkz( double rm, double dm, double amprms[21],
                      double $*$ra, double $*$da )
   }
   \sstarguments{
      \sstsubsection{
         rm = double (Given)
      }{
         Mean RA (radians).
      }
      \sstsubsection{
         dm = double (Given)
      }{
         Mean Dec (radians).
      }
      \sstsubsection{
         amprms = double[21] (Given)
      }{
         Star-independent mean-to-apparent parameters (see palMappa):
         (0-3)    not used
         (4-6)    heliocentric direction of the Earth (unit vector)
         (7)      not used
         (8-10)   abv: barycentric Earth velocity in units of c
         (11)     sqrt(1-v$\wedge$2) where v=modulus(abv)
         (12-20)  precession/nutation (3,3) matrix
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         Apparent RA (radians).
      }
      \sstsubsection{
         da = double $*$ (Returned)
      }{
         Apparent Dec (radians).
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The reference systems and timescales used are IAU 2006.

         \sstitem
          The mean place rm, dm and the vectors amprms[1-3] and amprms[4-6]
           are referred to the mean equinox and equator of the epoch
           specified when generating the precession/nutation matrix
           amprms[12-20].  In the call to palMappa (q.v.) normally used
           to populate amprms, this epoch is the first argument (eq).

         \sstitem
          The vector amprms(4-6) is referred to the mean equinox and
           equator of epoch eq.

         \sstitem
          Strictly speaking, the routine is not valid for solar-system
           sources, though the error will usually be extremely small.
           However, to prevent gross errors in the case where the
           position of the Sun is specified, the gravitational
           deflection term is restrained within about 920 arcsec of the
           centre of the Sun\texttt{'} s disc.  The term has a maximum value of
           about 1.85 arcsec at this radius, and decreases to zero as
           the centre of the disc is approached.
      }
   }
}
\sstroutine{
   palNut
}{
   Form the matrix of nutation
}{
   \sstdescription{
      Form the matrix of nutation for a given date using
      the IAU 2006 nutation model and palDeuler.
   }
   \sstinvocation{
      void palNut( double date, double rmatn[3][3] );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT as modified Julian date (JD-2400000.5)
      }
      \sstsubsection{
         rmatn = double [3][3] (Returned)
      }{
         Nutation matrix in the sense v(true)=rmatn $*$ v(mean)
         where v(true) is the star vector relative to the
         true equator and equinox of date and v(mean) is the
         star vector relative to the mean equator and equinox
         of date.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraNut06a via palNutc

         \sstitem
          The distinction between TDB and TT is negligible. For all but
           the most critical applications UTC is adequate.
      }
   }
}
\sstroutine{
   palNutc
}{
   Calculate nutation longitude \& obliquoty components
}{
   \sstdescription{
      Calculates the longitude $*$ obliquity components and mean obliquity
      using the SOFA/ERFA library.
   }
   \sstinvocation{
      void palNutc( double date, double $*$ dpsi, double $*$deps, double $*$eps0 );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT as modified Julian date (JD-2400000.5)
      }
      \sstsubsection{
         dpsi = double $*$ (Returned)
      }{
         Nutation in longitude
      }
      \sstsubsection{
         deps = double $*$ (Returned)
      }{
         Nutation in obliquity
      }
      \sstsubsection{
         eps0 = double $*$ (Returned)
      }{
         Mean obliquity.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Calls eraObl06 and eraNut06a and therefore uses the IAU 206
           precession/nutation model.

         \sstitem
          Note the change from SLA/F regarding the date. TT is used
           rather than TDB.
      }
   }
}
\sstroutine{
   palOap
}{
   Observed to apparent place
}{
   \sstdescription{
      Observed to apparent place.
   }
   \sstinvocation{
      void palOap ( const char $*$type, double ob1, double ob2, double date,
                    double dut, double elongm, double phim, double hm,
                    double xp, double yp, double tdk, double pmb,
                    double rh, double wl, double tlr,
                    double $*$rap, double $*$dap );
   }
   \sstarguments{
      \sstsubsection{
         type = const char $*$ (Given)
      }{
         Type of coordinates - \texttt{'} R\texttt{'} , \texttt{'} H\texttt{'}  or \texttt{'} A\texttt{'}  (see below)
      }
      \sstsubsection{
         ob1 = double (Given)
      }{
         Observed Az, HA or RA (radians; Az is N=0;E=90)
      }
      \sstsubsection{
         ob2 = double (Given)
      }{
         Observed ZD or Dec (radians)
      }
      \sstsubsection{
         date = double (Given)
      }{
         UTC date/time (Modified Julian Date, JD-2400000.5)
      }
      \sstsubsection{
         dut = double (Given)
      }{
         delta UT: UT1-UTC (UTC seconds)
      }
      \sstsubsection{
         elongm = double (Given)
      }{
         Mean longitude of the observer (radians, east $+$ve)
      }
      \sstsubsection{
         phim = double (Given)
      }{
         Mean geodetic latitude of the observer (radians)
      }
      \sstsubsection{
         hm = double (Given)
      }{
         Observer\texttt{'} s height above sea level (metres)
      }
      \sstsubsection{
         xp = double (Given)
      }{
         Polar motion x-coordinates (radians)
      }
      \sstsubsection{
         yp = double (Given)
      }{
         Polar motion y-coordinates (radians)
      }
      \sstsubsection{
         tdk = double (Given)
      }{
         Local ambient temperature (K; std=273.15)
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         Local atmospheric pressure (mb; std=1013.25)
      }
      \sstsubsection{
         rh = double (Given)
      }{
         Local relative humidity (in the range 0.0-1.0)
      }
      \sstsubsection{
         wl = double (Given)
      }{
         Effective wavelength (micron, e.g. 0.55)
      }
      \sstsubsection{
         tlr = double (Given)
      }{
         Tropospheric laps rate (K/metre, e.g. 0.0065)
      }
      \sstsubsection{
         rap = double $*$ (Given)
      }{
         Geocentric apparent right ascension
      }
      \sstsubsection{
         dap = double $*$ (Given)
      }{
         Geocentric apparent declination
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Only the first character of the TYPE argument is significant.
         \texttt{'} R\texttt{'}  or \texttt{'} r\texttt{'}  indicates that OBS1 and OBS2 are the observed right
         ascension and declination;  \texttt{'} H\texttt{'}  or \texttt{'} h\texttt{'}  indicates that they are
         hour angle (west $+$ve) and declination;  anything else (\texttt{'} A\texttt{'}  or
         \texttt{'} a\texttt{'}  is recommended) indicates that OBS1 and OBS2 are azimuth
         (north zero, east 90 deg) and zenith distance.  (Zenith
         distance is used rather than elevation in order to reflect the
         fact that no allowance is made for depression of the horizon.)

         \sstitem
          The accuracy of the result is limited by the corrections for
         refraction.  Providing the meteorological parameters are
         known accurately and there are no gross local effects, the
         predicted apparent RA,Dec should be within about 0.1 arcsec
         for a zenith distance of less than 70 degrees.  Even at a
         topocentric zenith distance of 90 degrees, the accuracy in
         elevation should be better than 1 arcmin;  useful results
         are available for a further 3 degrees, beyond which the
         palRefro routine returns a fixed value of the refraction.
         The complementary routines palAop (or palAopqk) and palOap
         (or palOapqk) are self-consistent to better than 1 micro-
         arcsecond all over the celestial sphere.

         \sstitem
          It is advisable to take great care with units, as even
         unlikely values of the input parameters are accepted and
         processed in accordance with the models used.

         \sstitem
          \texttt{"} Observed\texttt{"}  Az,El means the position that would be seen by a
         perfect theodolite located at the observer.  This is
         related to the observed HA,Dec via the standard rotation, using
         the geodetic latitude (corrected for polar motion), while the
         observed HA and RA are related simply through the local
         apparent ST.  \texttt{"} Observed\texttt{"}  RA,Dec or HA,Dec thus means the
         position that would be seen by a perfect equatorial located
         at the observer and with its polar axis aligned to the
         Earth\texttt{'} s axis of rotation (n.b. not to the refracted pole).
         By removing from the observed place the effects of
         atmospheric refraction and diurnal aberration, the
         geocentric apparent RA,Dec is obtained.

         \sstitem
          Frequently, mean rather than apparent RA,Dec will be required,
         in which case further transformations will be necessary.  The
         palAmp etc routines will convert the apparent RA,Dec produced
         by the present routine into an \texttt{"} FK5\texttt{"}  (J2000) mean place, by
         allowing for the Sun\texttt{'} s gravitational lens effect, annual
         aberration, nutation and precession.  Should \texttt{"} FK4\texttt{"}  (1950)
         coordinates be needed, the routines palFk524 etc will also
         need to be applied.

         \sstitem
          To convert to apparent RA,Dec the coordinates read from a
         real telescope, corrections would have to be applied for
         encoder zero points, gear and encoder errors, tube flexure,
         the position of the rotator axis and the pointing axis
         relative to it, non-perpendicularity between the mounting
         axes, and finally for the tilt of the azimuth or polar axis
         of the mounting (with appropriate corrections for mount
         flexures).  Some telescopes would, of course, exhibit other
         properties which would need to be accounted for at the
         appropriate point in the sequence.

         \sstitem
          This routine takes time to execute, due mainly to the rigorous
         integration used to evaluate the refraction.  For processing
         multiple stars for one location and time, call palAoppa once
         followed by one call per star to palOapqk.  Where a range of
         times within a limited period of a few hours is involved, and the
         highest precision is not required, call palAoppa once, followed
         by a call to palAoppat each time the time changes, followed by
         one call per star to palOapqk.

         \sstitem
          The DATE argument is UTC expressed as an MJD.  This is, strictly
         speaking, wrong, because of leap seconds.  However, as long as
         the delta UT and the UTC are consistent there are no
         difficulties, except during a leap second.  In this case, the
         start of the 61st second of the final minute should begin a new
         MJD day and the old pre-leap delta UT should continue to be used.
         As the 61st second completes, the MJD should revert to the start
         of the day as, simultaneously, the delta UTC changes by one
         second to its post-leap new value.

         \sstitem
          The delta UT (UT1-UTC) is tabulated in IERS circulars and
         elsewhere.  It increases by exactly one second at the end of
         each UTC leap second, introduced in order to keep delta UT
         within $+$/- 0.9 seconds.

         \sstitem
          IMPORTANT -- TAKE CARE WITH THE LONGITUDE SIGN CONVENTION.
         The longitude required by the present routine is east-positive,
         in accordance with geographical convention (and right-handed).
         In particular, note that the longitudes returned by the
         palOBS routine are west-positive, following astronomical
         usage, and must be reversed in sign before use in the present
         routine.

         \sstitem
          The polar coordinates XP,YP can be obtained from IERS
         circulars and equivalent publications.  The maximum amplitude
         is about 0.3 arcseconds.  If XP,YP values are unavailable,
         use XP=YP=0D0.  See page B60 of the 1988 Astronomical Almanac
         for a definition of the two angles.

         \sstitem
          The height above sea level of the observing station, HM,
         can be obtained from the Astronomical Almanac (Section J
         in the 1988 edition), or via the routine palOBS.  If P,
         the pressure in millibars, is available, an adequate
         estimate of HM can be obtained from the expression

      }
             HM $\sim$ -29.3$*$TSL$*$LOG(P/1013.25).

      where TSL is the approximate sea-level air temperature in K
      (see Astrophysical Quantities, C.W.Allen, 3rd edition,
      section 52).  Similarly, if the pressure P is not known,
      it can be estimated from the height of the observing
      station, HM, as follows:

             P $\sim$ 1013.25$*$EXP(-HM/(29.3$*$TSL)).

      Note, however, that the refraction is nearly proportional to the
      pressure and that an accurate P value is important for precise
      work.

      \sstitemlist{

         \sstitem
          The azimuths etc. used by the present routine are with respect
         to the celestial pole.  Corrections from the terrestrial pole
         can be computed using palPolmo.
      }
   }
}
\sstroutine{
   palOapqk
}{
   Quick observed to apparent place
}{
   \sstdescription{
      type = const char $*$ (Given)
         Type of coordinates - \texttt{'} R\texttt{'} , \texttt{'} H\texttt{'}  or \texttt{'} A\texttt{'}  (see below)
      ob1 = double (Given)
         Observed Az, HA or RA (radians; Az is N=0;E=90)
      ob2 = double (Given)
         Observed ZD or Dec (radians)
      aoprms = const double [14] (Given)
         Star-independent apparent-to-observed parameters.
         See palAopqk for details.
      rap = double $*$ (Given)
         Geocentric apparent right ascension
      dap = double $*$ (Given)
         Geocentric apparent declination
   }
   \sstinvocation{
      void palOapqk ( const char $*$type, double ob1, double ob2,
                      const  double aoprms[14], double $*$rap, double $*$dap );
   }
   \sstarguments{
      \sstsubsection{
         Quick observed to apparent place.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Only the first character of the TYPE argument is significant.
         \texttt{'} R\texttt{'}  or \texttt{'} r\texttt{'}  indicates that OBS1 and OBS2 are the observed right
         ascension and declination;  \texttt{'} H\texttt{'}  or \texttt{'} h\texttt{'}  indicates that they are
         hour angle (west $+$ve) and declination;  anything else (\texttt{'} A\texttt{'}  or
         \texttt{'} a\texttt{'}  is recommended) indicates that OBS1 and OBS2 are azimuth
         (north zero, east 90 deg) and zenith distance.  (Zenith distance
         is used rather than elevation in order to reflect the fact that
         no allowance is made for depression of the horizon.)

         \sstitem
          The accuracy of the result is limited by the corrections for
         refraction.  Providing the meteorological parameters are
         known accurately and there are no gross local effects, the
         predicted apparent RA,Dec should be within about 0.1 arcsec
         for a zenith distance of less than 70 degrees.  Even at a
         topocentric zenith distance of 90 degrees, the accuracy in
         elevation should be better than 1 arcmin;  useful results
         are available for a further 3 degrees, beyond which the
         palREFRO routine returns a fixed value of the refraction.
         The complementary routines palAop (or palAopqk) and palOap
         (or palOapqk) are self-consistent to better than 1 micro-
         arcsecond all over the celestial sphere.

         \sstitem
          It is advisable to take great care with units, as even
         unlikely values of the input parameters are accepted and
         processed in accordance with the models used.

         \sstitem
          \texttt{"} Observed\texttt{"}  Az,El means the position that would be seen by a
         perfect theodolite located at the observer.  This is
         related to the observed HA,Dec via the standard rotation, using
         the geodetic latitude (corrected for polar motion), while the
         observed HA and RA are related simply through the local
         apparent ST.  \texttt{"} Observed\texttt{"}  RA,Dec or HA,Dec thus means the
         position that would be seen by a perfect equatorial located
         at the observer and with its polar axis aligned to the
         Earth\texttt{'} s axis of rotation (n.b. not to the refracted pole).
         By removing from the observed place the effects of
         atmospheric refraction and diurnal aberration, the
         geocentric apparent RA,Dec is obtained.

         \sstitem
          Frequently, mean rather than apparent RA,Dec will be required,
         in which case further transformations will be necessary.  The
         palAmp etc routines will convert the apparent RA,Dec produced
         by the present routine into an \texttt{"} FK5\texttt{"}  (J2000) mean place, by
         allowing for the Sun\texttt{'} s gravitational lens effect, annual
         aberration, nutation and precession.  Should \texttt{"} FK4\texttt{"}  (1950)
         coordinates be needed, the routines palFk524 etc will also
         need to be applied.

         \sstitem
          To convert to apparent RA,Dec the coordinates read from a
         real telescope, corrections would have to be applied for
         encoder zero points, gear and encoder errors, tube flexure,
         the position of the rotator axis and the pointing axis
         relative to it, non-perpendicularity between the mounting
         axes, and finally for the tilt of the azimuth or polar axis
         of the mounting (with appropriate corrections for mount
         flexures).  Some telescopes would, of course, exhibit other
         properties which would need to be accounted for at the
         appropriate point in the sequence.

         \sstitem
          The star-independent apparent-to-observed-place parameters
         in AOPRMS may be computed by means of the palAoppa routine.
         If nothing has changed significantly except the time, the
         palAoppat routine may be used to perform the requisite
         partial recomputation of AOPRMS.

         \sstitem
          The azimuths etc used by the present routine are with respect
         to the celestial pole.  Corrections from the terrestrial pole
         can be computed using palPolmo.
      }
   }
}
\sstroutine{
   palObs
}{
   Parameters of selected ground-based observing stations
}{
   \sstdescription{
      Station numbers, identifiers, names and other details are
      subject to change and should not be hardwired into
      application programs.

      All characters in \texttt{"} c\texttt{"}  up to the first space are
      checked;  thus an abbreviated ID will return the parameters
      for the first station in the list which matches the
      abbreviation supplied, and no station in the list will ever
      contain embedded spaces. \texttt{"} c\texttt{"}  must not have leading spaces.

      IMPORTANT -- BEWARE OF THE LONGITUDE SIGN CONVENTION.  The
      longitude returned by palOBS (and SLA\_OBS) is west-positive in accordance
      with astronomical usage.  However, this sign convention is
      left-handed and is the opposite of the one used by geographers;
      elsewhere in PAL the preferable east-positive convention is
      used.  In particular, note that for use in palAop, palAoppa
      and palOap the sign of the longitude must be reversed.

      Users are urged to inform the author of any improvements
      they would like to see made.  For example:

          typographical corrections
          more accurate parameters
          better station identifiers or names
          additional stations
   }
   \sstinvocation{
      int palObs( size\_t n, const char $*$ c,
                  char $*$ ident, size\_t identlen,
                  char $*$ name, size\_t namelen,
                  double $*$ w, double $*$ p, double $*$ h );
   }
   \sstarguments{
      \sstsubsection{
         n = size\_t (Given)
      }{
         Number specifying the observing station. If 0
         the identifier in \texttt{"} c\texttt{"}  is used to determine the
         observing station to use.
      }
      \sstsubsection{
         c = const char $*$ (Given)
      }{
         Identifier specifying the observing station for
         which the parameters should be returned. Only used
         if n is 0. Can be NULL for n$>$0. Case insensitive.
      }
      \sstsubsection{
         ident = char $*$ (Returned)
      }{
         Identifier of the observing station selected. Will be
         identical to \texttt{"} c\texttt{"}  if n==0. Unchanged if \texttt{"} n\texttt{"}  or \texttt{"} c\texttt{"}
         do not match an observing station. Should be at least
         11 characters (including the trailing nul).
      }
      \sstsubsection{
         identlen = size\_t (Given)
      }{
         Size of the buffer \texttt{"} ident\texttt{"}  including trailing nul.
      }
      \sstsubsection{
         name = char $*$ (Returned)
      }{
         Full name of the specified observing station. Contains \texttt{"} ?\texttt{"}
         if \texttt{"} n\texttt{"}  or \texttt{"} c\texttt{"}  did not correspond to a valid station. Should
         be at least 41 characters (including the trailing nul).
      }
      \sstsubsection{
         w = double $*$ (Returned)
      }{
         Longitude (radians, West $+$ve). Unchanged if observing
         station could not be identified.
      }
      \sstsubsection{
         p = double $*$ (Returned)
      }{
         Geodetic latitude (radians, North $+$ve). Unchanged if observing
         station could not be identified.
      }
      \sstsubsection{
         h = double $*$ (Returned)
      }{
         Height above sea level (metres). Unchanged if observing
         station could not be identified.
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         palObs = int
      }{
         0 if an observing station was returned. -1 if no match was
         found.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Differs from the SLA interface in that the output short name
           is not the same variable as the input short name. This simplifies
           consting. Additionally the size of the output buffers are now
           specified in the API and a status integer is returned.
      }
   }
}
\sstroutine{
   palPa
}{
   HA, Dec to Parallactic Angle
}{
   \sstdescription{
      Converts HA, Dec to Parallactic Angle.
   }
   \sstinvocation{
      double palPa( double ha, double dec, double phi );
   }
   \sstarguments{
      \sstsubsection{
         ha = double (Given)
      }{
         Hour angle in radians (Geocentric apparent)
      }
      \sstsubsection{
         dec = double (Given)
      }{
         Declination in radians (Geocentric apparent)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observatory latitude in radians (geodetic)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         palPa = double
      }{
         Parallactic angle in the range -pi to $+$pi.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The parallactic angle at a point in the sky is the position
           angle of the vertical, i.e. the angle between the direction to
           the pole and to the zenith.  In precise applications care must
           be taken only to use geocentric apparent HA,Dec and to consider
           separately the effects of atmospheric refraction and telescope
           mount errors.

         \sstitem
          At the pole a zero result is returned.
      }
   }
}
\sstroutine{
   palPcd
}{
   Apply pincushion/barrel distortion to a tangent-plane [x,y]
}{
   \sstdescription{
      Applies pincushion and barrel distortion to a tangent
      plane coordinate.
   }
   \sstinvocation{
      palPcd( double disco, double $*$ x, double $*$ y );
   }
   \sstarguments{
      \sstsubsection{
         disco = double (Given)
      }{
         Pincushion/barrel distortion coefficient.
      }
      \sstsubsection{
         x = double $*$ (Given \& Returned)
      }{
         On input the tangent-plane X coordinate, on output
         the distorted X coordinate.
      }
      \sstsubsection{
         y = double $*$ (Given \& Returned)
      }{
         On input the tangent-plane Y coordinate, on output
         the distorted Y coordinate.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The distortion is of the form RP = R$*$(1 $+$ C$*$R$*$$*$2), where R is
           the radial distance from the tangent point, C is the DISCO
           argument, and RP is the radial distance in the presence of
           the distortion.

         \sstitem
          For pincushion distortion, C is $+$ve;  for barrel distortion,
           C is -ve.

         \sstitem
          For X,Y in units of one projection radius (in the case of
           a photographic plate, the focal length), the following
           DISCO values apply:

      }
            Geometry          DISCO

            astrograph         0.0
            Schmidt           -0.3333
            AAT PF doublet  $+$147.069
            AAT PF triplet  $+$178.585
            AAT f/8          $+$21.20
            JKT f/8          $+$13.32
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          There is a companion routine, palUnpcd, which performs the
           inverse operation.
      }
   }
}
\sstroutine{
   palPertel
}{
   Update elements by applying planetary perturbations
}{
   \sstdescription{
      Update the osculating orbital elements of an asteroid or comet by
      applying planetary perturbations.
   }
   \sstinvocation{
      void palPertel (int jform, double date0, double date1,
                      double epoch0, double orbi0, double anode0,
                      double perih0, double aorq0, double e0, double am0,
                      double $*$epoch1, double $*$orbi1, double $*$anode1,
                      double $*$perih1, double $*$aorq1, double $*$e1, double $*$am1,
                      int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         jform = int (Given)
      }{
         Element set actually returned (1-3; Note 6)
      }
      \sstsubsection{
         date0 = double (Given)
      }{
         Date of osculation (TT MJD) for the given elements.
      }
      \sstsubsection{
         date1 = double (Given)
      }{
         Date of osculation (TT MJD) for the updated elements.
      }
      \sstsubsection{
         epoch0 = double (Given)
      }{
         Epoch of elements (TT MJD)
      }
      \sstsubsection{
         orbi0 = double (Given)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode0 = double (Given)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih0 = double (Given)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq0 = double (Given)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e0 = double (Given)
      }{
         eccentricity
      }
      \sstsubsection{
         am0 = double (Given)
      }{
         mean anomaly (radians, JFORM=2 only)
      }
      \sstsubsection{
         epoch1 = double $*$ (Returned)
      }{
         Epoch of elements (TT MJD)
      }
      \sstsubsection{
         orbi1 = double $*$ (Returned)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode1 = double $*$ (Returned)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih1 = double $*$ (Returned)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq1 = double $*$ (Returned)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e1 = double $*$ (Returned)
      }{
         eccentricity
      }
      \sstsubsection{
         am1 = double $*$ (Returned)
      }{
         mean anomaly (radians, JFORM=2 only)
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:
         \sstitemlist{

            \sstitem
               $+$102 = warning, distant epoch

            \sstitem
               $+$101 = warning, large timespan ( $>$ 100 years)

            \sstitem
               $+$1 to $+$10 = coincident with planet (Note 6)

            \sstitem
               0 = OK

            \sstitem
               -1 = illegal JFORM

            \sstitem
               -2 = illegal E0

            \sstitem
               -3 = illegal AORQ0

            \sstitem
               -4 = internal error

            \sstitem
               -5 = numerical error
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Two different element-format options are available:

      }
        Option JFORM=2, suitable for minor planets:

        EPOCH   = epoch of elements (TT MJD)
        ORBI    = inclination i (radians)
        ANODE   = longitude of the ascending node, big omega (radians)
        PERIH   = argument of perihelion, little omega (radians)
        AORQ    = mean distance, a (AU)
        E       = eccentricity, e
        AM      = mean anomaly M (radians)

        Option JFORM=3, suitable for comets:

        EPOCH   = epoch of perihelion (TT MJD)
        ORBI    = inclination i (radians)
        ANODE   = longitude of the ascending node, big omega (radians)
        PERIH   = argument of perihelion, little omega (radians)
        AORQ    = perihelion distance, q (AU)
        E       = eccentricity, e

      \sstitemlist{

         \sstitem
          DATE0, DATE1, EPOCH0 and EPOCH1 are all instants of time in
           the TT timescale (formerly Ephemeris Time, ET), expressed
           as Modified Julian Dates (JD-2400000.5).

      }
        DATE0 is the instant at which the given (i.e. unperturbed)
        osculating elements are correct.

        DATE1 is the specified instant at which the updated osculating
        elements are correct.

        EPOCH0 and EPOCH1 will be the same as DATE0 and DATE1
        (respectively) for the JFORM=2 case, normally used for minor
        planets.  For the JFORM=3 case, the two epochs will refer to
        perihelion passage and so will not, in general, be the same as
        DATE0 and/or DATE1 though they may be similar to one another.
      \sstitemlist{

         \sstitem
          The elements are with respect to the J2000 ecliptic and equinox.

         \sstitem
          Unused elements (AM0 and AM1 for JFORM=3) are not accessed.

         \sstitem
          See the palPertue routine for details of the algorithm used.

         \sstitem
          This routine is not intended to be used for major planets, which
           is why JFORM=1 is not available and why there is no opportunity
           to specify either the longitude of perihelion or the daily
           motion.  However, if JFORM=2 elements are somehow obtained for a
           major planet and supplied to the routine, sensible results will,
           in fact, be produced.  This happens because the palPertue routine
           that is called to perform the calculations checks the separation
           between the body and each of the planets and interprets a
           suspiciously small value (0.001 AU) as an attempt to apply it to
           the planet concerned.  If this condition is detected, the
           contribution from that planet is ignored, and the status is set to
           the planet number (1-10 = Mercury, Venus, EMB, Mars, Jupiter,
           Saturn, Uranus, Neptune, Earth, Moon) as a warning.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          Sterne, Theodore E., \texttt{"} An Introduction to Celestial Mechanics\texttt{"} ,
           Interscience Publishers Inc., 1960.  Section 6.7, p199.
      }
   }
}
\sstroutine{
   palPertue
}{
   Update the universal elements by applying planetary perturbations
}{
   \sstdescription{
      Update the universal elements of an asteroid or comet by applying
      planetary perturbations.
   }
   \sstinvocation{
      void palPertue( double date, double u[13], int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         Final epoch (TT MJD) for the update elements.
      }
      \sstsubsection{
         u = const double [13] (Given \& Returned)
      }{
         Universal orbital elements (Note 1)
             (0)  combined mass (M$+$m)
             (1)  total energy of the orbit (alpha)
             (2)  reference (osculating) epoch (t0)
           (3-5)  position at reference epoch (r0)
           (6-8)  velocity at reference epoch (v0)
             (9)  heliocentric distance at reference epoch
            (10)  r0.v0
            (11)  date (t)
            (12)  universal eccentric anomaly (psi) of date, approx
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:
                    $+$102 = warning, distant epoch
                    $+$101 = warning, large timespan ( $>$ 100 years)
               $+$1 to $+$10 = coincident with major planet (Note 5)
                       0 = OK
         \sstitemlist{

            \sstitem
                          1 = numerical error
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The \texttt{"} universal\texttt{"}  elements are those which define the orbit for the
           purposes of the method of universal variables (see reference 2).
           They consist of the combined mass of the two bodies, an epoch,
           and the position and velocity vectors (arbitrary reference frame)
           at that epoch.  The parameter set used here includes also various
           quantities that can, in fact, be derived from the other
           information.  This approach is taken to avoiding unnecessary
           computation and loss of accuracy.  The supplementary quantities
           are (i) alpha, which is proportional to the total energy of the
           orbit, (ii) the heliocentric distance at epoch, (iii) the
           outwards component of the velocity at the given epoch, (iv) an
           estimate of psi, the \texttt{"} universal eccentric anomaly\texttt{"}  at a given
           date and (v) that date.

         \sstitem
          The universal elements are with respect to the J2000 equator and
           equinox.

         \sstitem
          The epochs DATE, U(3) and U(12) are all Modified Julian Dates
           (JD-2400000.5).

         \sstitem
          The algorithm is a simplified form of Encke\texttt{'} s method.  It takes as
           a basis the unperturbed motion of the body, and numerically
           integrates the perturbing accelerations from the major planets.
           The expression used is essentially Sterne\texttt{'} s 6.7-2 (reference 1).
           Everhart and Pitkin (reference 2) suggest rectifying the orbit at
           each integration step by propagating the new perturbed position
           and velocity as the new universal variables.  In the present
           routine the orbit is rectified less frequently than this, in order
           to gain a slight speed advantage.  However, the rectification is
           done directly in terms of position and velocity, as suggested by
           Everhart and Pitkin, bypassing the use of conventional orbital
           elements.

      }
        The f(q) part of the full Encke method is not used.  The purpose
        of this part is to avoid subtracting two nearly equal quantities
        when calculating the \texttt{"} indirect member\texttt{"} , which takes account of the
        small change in the Sun\texttt{'} s attraction due to the slightly displaced
        position of the perturbed body.  A simpler, direct calculation in
        double precision proves to be faster and not significantly less
        accurate.

        Apart from employing a variable timestep, and occasionally
        \texttt{"} rectifying the orbit\texttt{"}  to keep the indirect member small, the
        integration is done in a fairly straightforward way.  The
        acceleration estimated for the middle of the timestep is assumed
        to apply throughout that timestep;  it is also used in the
        extrapolation of the perturbations to the middle of the next
        timestep, to predict the new disturbed position.  There is no
        iteration within a timestep.

        Measures are taken to reach a compromise between execution time
        and accuracy.  The starting-point is the goal of achieving
        arcsecond accuracy for ordinary minor planets over a ten-year
        timespan.  This goal dictates how large the timesteps can be,
        which in turn dictates how frequently the unperturbed motion has
        to be recalculated from the osculating elements.

        Within predetermined limits, the timestep for the numerical
        integration is varied in length in inverse proportion to the
        magnitude of the net acceleration on the body from the major
        planets.

        The numerical integration requires estimates of the major-planet
        motions.  Approximate positions for the major planets (Pluto
        alone is omitted) are obtained from the routine palPlanet.  Two
        levels of interpolation are used, to enhance speed without
        significantly degrading accuracy.  At a low frequency, the routine
        palPlanet is called to generate updated position$+$velocity \texttt{"} state
        vectors\texttt{"} .  The only task remaining to be carried out at the full
        frequency (i.e. at each integration step) is to use the state
        vectors to extrapolate the planetary positions.  In place of a
        strictly linear extrapolation, some allowance is made for the
        curvature of the orbit by scaling back the radius vector as the
        linear extrapolation goes off at a tangent.

        Various other approximations are made.  For example, perturbations
        by Pluto and the minor planets are neglected and relativistic
        effects are not taken into account.

        In the interests of simplicity, the background calculations for
        the major planets are carried out en masse.  The mean elements and
        state vectors for all the planets are refreshed at the same time,
        without regard for orbit curvature, mass or proximity.

        The Earth-Moon system is treated as a single body when the body is
        distant but as separate bodies when closer to the EMB than the
        parameter RNE, which incurs a time penalty but improves accuracy
        for near-Earth objects.

      \sstitemlist{

         \sstitem
          This routine is not intended to be used for major planets.
           However, if major-planet elements are supplied, sensible results
           will, in fact, be produced.  This happens because the routine
           checks the separation between the body and each of the planets and
           interprets a suspiciously small value (0.001 AU) as an attempt to
           apply the routine to the planet concerned.  If this condition is
           detected, the contribution from that planet is ignored, and the
           status is set to the planet number (1-10 = Mercury, Venus, EMB,
           Mars, Jupiter, Saturn, Uranus, Neptune, Earth, Moon) as a warning.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          Sterne, Theodore E., \texttt{"} An Introduction to Celestial Mechanics\texttt{"} ,
           Interscience Publishers Inc., 1960.  Section 6.7, p199.

         \sstitem
          Everhart, E. \& Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
      }
   }
}
\sstroutine{
   palPlanel
}{
   Transform conventional elements into position and velocity
}{
   \sstdescription{
      Heliocentric position and velocity of a planet, asteroid or comet,
      starting from orbital elements.
   }
   \sstinvocation{
      void palPlanel ( double date, int jform, double epoch, double orbinc,
                       double anode, double perih, double aorq, double e,
                       double aorl, double dm, double pv[6], int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         Epoch (TT MJD) of osculation (Note 1)
      }
      \sstsubsection{
         jform = int (Given)
      }{
         Element set actually returned (1-3; Note 3)
      }
      \sstsubsection{
         epoch = double (Given)
      }{
         Epoch of elements (TT MJD) (Note 4)
      }
      \sstsubsection{
         orbinc = double (Given)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode = double (Given)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih = double (Given)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq = double (Given)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e = double (Given)
      }{
         eccentricity
      }
      \sstsubsection{
         aorl = double (Given)
      }{
         mean anomaly or longitude (radians, JFORM=1,2 only)
      }
      \sstsubsection{
         dm = double (Given)
      }{
         daily motion (radians, JFORM=1 only)
      }
      \sstsubsection{
         u = double [13] (Returned)
      }{
         Universal orbital elements (Note 1)
             (0)  combined mass (M$+$m)
             (1)  total energy of the orbit (alpha)
             (2)  reference (osculating) epoch (t0)
           (3-5)  position at reference epoch (r0)
           (6-8)  velocity at reference epoch (v0)
             (9)  heliocentric distance at reference epoch
            (10)  r0.v0
            (11)  date (t)
            (12)  universal eccentric anomaly (psi) of date, approx
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:  0 = OK
         \sstitemlist{

            \sstitem
                   -1 = illegal JFORM

            \sstitem
                   -2 = illegal E

            \sstitem
                   -3 = illegal AORQ

            \sstitem
                   -4 = illegal DM

            \sstitem
                   -5 = numerical error
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          DATE is the instant for which the prediction is required.  It is
           in the TT timescale (formerly Ephemeris Time, ET) and is a
           Modified Julian Date (JD-2400000.5).

         \sstitem
          The elements are with respect to the J2000 ecliptic and equinox.

         \sstitem
          A choice of three different element-set options is available:

      }
        Option JFORM = 1, suitable for the major planets:

          EPOCH  = epoch of elements (TT MJD)
          ORBINC = inclination i (radians)
          ANODE  = longitude of the ascending node, big omega (radians)
          PERIH  = longitude of perihelion, curly pi (radians)
          AORQ   = mean distance, a (AU)
          E      = eccentricity, e (range 0 to $<$1)
          AORL   = mean longitude L (radians)
          DM     = daily motion (radians)

        Option JFORM = 2, suitable for minor planets:

          EPOCH  = epoch of elements (TT MJD)
          ORBINC = inclination i (radians)
          ANODE  = longitude of the ascending node, big omega (radians)
          PERIH  = argument of perihelion, little omega (radians)
          AORQ   = mean distance, a (AU)
          E      = eccentricity, e (range 0 to $<$1)
          AORL   = mean anomaly M (radians)

        Option JFORM = 3, suitable for comets:

          EPOCH  = epoch of elements and perihelion (TT MJD)
          ORBINC = inclination i (radians)
          ANODE  = longitude of the ascending node, big omega (radians)
          PERIH  = argument of perihelion, little omega (radians)
          AORQ   = perihelion distance, q (AU)
          E      = eccentricity, e (range 0 to 10)

        Unused arguments (DM for JFORM=2, AORL and DM for JFORM=3) are not
        accessed.
      \sstitemlist{

         \sstitem
          Each of the three element sets defines an unperturbed heliocentric
           orbit.  For a given epoch of observation, the position of the body
           in its orbit can be predicted from these elements, which are
           called \texttt{"} osculating elements\texttt{"} , using standard two-body analytical
           solutions.  However, due to planetary perturbations, a given set
           of osculating elements remains usable for only as long as the
           unperturbed orbit that it describes is an adequate approximation
           to reality.  Attached to such a set of elements is a date called
           the \texttt{"} osculating epoch\texttt{"} , at which the elements are, momentarily,
           a perfect representation of the instantaneous position and
           velocity of the body.

      }
        Therefore, for any given problem there are up to three different
        epochs in play, and it is vital to distinguish clearly between
        them:

        . The epoch of observation:  the moment in time for which the
          position of the body is to be predicted.

        . The epoch defining the position of the body:  the moment in time
          at which, in the absence of purturbations, the specified
          position (mean longitude, mean anomaly, or perihelion) is
          reached.

        . The osculating epoch:  the moment in time at which the given
          elements are correct.

        For the major-planet and minor-planet cases it is usual to make
        the epoch that defines the position of the body the same as the
        epoch of osculation.  Thus, only two different epochs are
        involved:  the epoch of the elements and the epoch of observation.

        For comets, the epoch of perihelion fixes the position in the
        orbit and in general a different epoch of osculation will be
        chosen.  Thus, all three types of epoch are involved.

        For the present routine:

        . The epoch of observation is the argument DATE.

        . The epoch defining the position of the body is the argument
          EPOCH.

        . The osculating epoch is not used and is assumed to be close
          enough to the epoch of observation to deliver adequate accuracy.
          If not, a preliminary call to palPertel may be used to update
          the element-set (and its associated osculating epoch) by
          applying planetary perturbations.
      \sstitemlist{

         \sstitem
          The reference frame for the result is with respect to the mean
           equator and equinox of epoch J2000.

         \sstitem
          The algorithm was originally adapted from the EPHSLA program of
           D.H.P.Jones (private communication, 1996).  The method is based
           on Stumpff\texttt{'} s Universal Variables.
      }
   }
   \sstdiytopic{
      See Also
   }{
      Everhart, E. \& Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
   }
}
\sstroutine{
   palPlanet
}{
   Approximate heliocentric position and velocity of major planet
}{
   \sstdescription{
      Calculates the approximate heliocentric position and velocity of
      the specified major planet.
   }
   \sstinvocation{
      void palPlanet ( double date, int np, double pv[6], int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TDB Modified Julian Date (JD-2400000.5).
      }
      \sstsubsection{
         np = int (Given)
      }{
         planet (1=Mercury, 2=Venus, 3=EMB, 4=Mars,
                 5=Jupiter, 6=Saturn, 7=Uranus, 8=Neptune)
      }
      \sstsubsection{
         pv = double [6] (Returned)
      }{
         heliocentric x,y,z,xdot,ydot,zdot, J2000, equatorial triad
         in units AU and AU/s.
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         \sstitemlist{

            \sstitem
             -2 = solution didn\texttt{'} t converge.

            \sstitem
             -1 = illegal np (1-8)

            \sstitem
             0 = OK

            \sstitem
             $+$1 = warning: year outside 1000-3000
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          See SOFA/ERFA eraPlan94 for details

         \sstitem
          Note that Pluto is supported in SLA/F but not in this routine

         \sstitem
          Status -2 is equivalent to eraPlan94 status $+$2.

         \sstitem
          Note that velocity units here match the SLA/F documentation.
      }
   }
}
\sstroutine{
   palPlante
}{
   Topocentric RA,Dec of a Solar-System object from heliocentric orbital elements
}{
   \sstdescription{
      Topocentric apparent RA,Dec of a Solar-System object whose
      heliocentric orbital elements are known.
   }
   \sstinvocation{
      void palPlante ( double date, double elong, double phi, int jform,
                       double epoch, double orbinc, double anode, double perih,
                       double aorq, double e, double aorl, double dm,
                       double $*$ra, double $*$dec, double $*$r, int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT MJD of observation (JD-2400000.5)
      }
      \sstsubsection{
         elong = double (Given)
      }{
         Observer\texttt{'} s east longitude (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observer\texttt{'} s geodetic latitude (radians)
      }
      \sstsubsection{
         jform = int (Given)
      }{
         Element set actually returned (1-3; Note 6)
      }
      \sstsubsection{
         epoch = double (Given)
      }{
         Epoch of elements (TT MJD)
      }
      \sstsubsection{
         orbinc = double (Given)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode = double (Given)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih = double (Given)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq = double (Given)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e = double (Given)
      }{
         eccentricity
      }
      \sstsubsection{
         aorl = double (Given)
      }{
         mean anomaly or longitude (radians, JFORM=1,2 only)
      }
      \sstsubsection{
         dm = double (Given)
      }{
         daily motion (radians, JFORM=1 only)
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         Topocentric apparent RA (radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Topocentric apparent Dec (radians)
      }
      \sstsubsection{
         r = double $*$ (Returned)
      }{
         Distance from observer (AU)
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status: 0 = OK
         \sstitemlist{

            \sstitem
                  -1 = illegal jform

            \sstitem
                  -2 = illegal e

            \sstitem
                  -3 = illegal aorq

            \sstitem
                  -4 = illegal dm

            \sstitem
                  -5 = numerical error
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          DATE is the instant for which the prediction is required.  It is
           in the TT timescale (formerly Ephemeris Time, ET) and is a
           Modified Julian Date (JD-2400000.5).

         \sstitem
          The longitude and latitude allow correction for geocentric
           parallax.  This is usually a small effect, but can become
           important for near-Earth asteroids.  Geocentric positions can be
           generated by appropriate use of routines palEpv (or palEvp) and
           palUe2pv.

         \sstitem
          The elements are with respect to the J2000 ecliptic and equinox.

         \sstitem
          A choice of three different element-set options is available:

      }
        Option JFORM = 1, suitable for the major planets:

          EPOCH  = epoch of elements (TT MJD)
          ORBINC = inclination i (radians)
          ANODE  = longitude of the ascending node, big omega (radians)
          PERIH  = longitude of perihelion, curly pi (radians)
          AORQ   = mean distance, a (AU)
          E      = eccentricity, e (range 0 to $<$1)
          AORL   = mean longitude L (radians)
          DM     = daily motion (radians)

        Option JFORM = 2, suitable for minor planets:

          EPOCH  = epoch of elements (TT MJD)
          ORBINC = inclination i (radians)
          ANODE  = longitude of the ascending node, big omega (radians)
          PERIH  = argument of perihelion, little omega (radians)
          AORQ   = mean distance, a (AU)
          E      = eccentricity, e (range 0 to $<$1)
          AORL   = mean anomaly M (radians)

        Option JFORM = 3, suitable for comets:

          EPOCH  = epoch of elements and perihelion (TT MJD)
          ORBINC = inclination i (radians)
          ANODE  = longitude of the ascending node, big omega (radians)
          PERIH  = argument of perihelion, little omega (radians)
          AORQ   = perihelion distance, q (AU)
          E      = eccentricity, e (range 0 to 10)

        Unused arguments (DM for JFORM=2, AORL and DM for JFORM=3) are not
        accessed.
      \sstitemlist{

         \sstitem
          Each of the three element sets defines an unperturbed heliocentric
           orbit.  For a given epoch of observation, the position of the body
           in its orbit can be predicted from these elements, which are
           called \texttt{"} osculating elements\texttt{"} , using standard two-body analytical
           solutions.  However, due to planetary perturbations, a given set
           of osculating elements remains usable for only as long as the
           unperturbed orbit that it describes is an adequate approximation
           to reality.  Attached to such a set of elements is a date called
           the \texttt{"} osculating epoch\texttt{"} , at which the elements are, momentarily,
           a perfect representation of the instantaneous position and
           velocity of the body.

      }
        Therefore, for any given problem there are up to three different
        epochs in play, and it is vital to distinguish clearly between
        them:

        . The epoch of observation:  the moment in time for which the
          position of the body is to be predicted.

        . The epoch defining the position of the body:  the moment in time
          at which, in the absence of purturbations, the specified
          position (mean longitude, mean anomaly, or perihelion) is
          reached.

        . The osculating epoch:  the moment in time at which the given
          elements are correct.

        For the major-planet and minor-planet cases it is usual to make
        the epoch that defines the position of the body the same as the
        epoch of osculation.  Thus, only two different epochs are
        involved:  the epoch of the elements and the epoch of observation.

        For comets, the epoch of perihelion fixes the position in the
        orbit and in general a different epoch of osculation will be
        chosen.  Thus, all three types of epoch are involved.

        For the present routine:

        . The epoch of observation is the argument DATE.

        . The epoch defining the position of the body is the argument
          EPOCH.

        . The osculating epoch is not used and is assumed to be close
          enough to the epoch of observation to deliver adequate accuracy.
          If not, a preliminary call to palPertel may be used to update
          the element-set (and its associated osculating epoch) by
          applying planetary perturbations.
      \sstitemlist{

         \sstitem
          Two important sources for orbital elements are Horizons, operated
           by the Jet Propulsion Laboratory, Pasadena, and the Minor Planet
           Center, operated by the Center for Astrophysics, Harvard.

      }
        The JPL Horizons elements (heliocentric, J2000 ecliptic and
        equinox) correspond to PAL/SLALIB arguments as follows.

         Major planets:

          JFORM  = 1
          EPOCH  = JDCT-2400000.5
          ORBINC = IN (in radians)
          ANODE  = OM (in radians)
          PERIH  = OM$+$W (in radians)
          AORQ   = A
          E      = EC
          AORL   = MA$+$OM$+$W (in radians)
          DM     = N (in radians)

          Epoch of osculation = JDCT-2400000.5

         Minor planets:

          JFORM  = 2
          EPOCH  = JDCT-2400000.5
          ORBINC = IN (in radians)
          ANODE  = OM (in radians)
          PERIH  = W (in radians)
          AORQ   = A
          E      = EC
          AORL   = MA (in radians)

          Epoch of osculation = JDCT-2400000.5

         Comets:

          JFORM  = 3
          EPOCH  = Tp-2400000.5
          ORBINC = IN (in radians)
          ANODE  = OM (in radians)
          PERIH  = W (in radians)
          AORQ   = QR
          E      = EC

          Epoch of osculation = JDCT-2400000.5

       The MPC elements correspond to SLALIB arguments as follows.

         Minor planets:

          JFORM  = 2
          EPOCH  = Epoch-2400000.5
          ORBINC = Incl. (in radians)
          ANODE  = Node (in radians)
          PERIH  = Perih. (in radians)
          AORQ   = a
          E      = e
          AORL   = M (in radians)

          Epoch of osculation = Epoch-2400000.5

        Comets:

          JFORM  = 3
          EPOCH  = T-2400000.5
          ORBINC = Incl. (in radians)
          ANODE  = Node. (in radians)
          PERIH  = Perih. (in radians)
          AORQ   = q
          E      = e

          Epoch of osculation = Epoch-2400000.5
   }
}
\sstroutine{
   palPlantu
}{
   Topocentric RA,Dec of a Solar-System object from universal elements
}{
   \sstdescription{
      Topocentric apparent RA,Dec of a Solar-System object whose
      heliocentric universal elements are known.
   }
   \sstinvocation{
      void palPlantu ( double date, double elong, double phi, const double u[13],
                       double $*$ra, double $*$dec, double $*$r, int $*$jstat ) \{
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT MJD of observation (JD-2400000.5)
      }
      \sstsubsection{
         elong = double (Given)
      }{
         Observer\texttt{'} s east longitude (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observer\texttt{'} s geodetic latitude (radians)
      }
      \sstsubsection{
         u = const double [13] (Given)
      }{
         Universal orbital elements
         \sstitemlist{

            \sstitem
               (0)  combined mass (M$+$m)

            \sstitem
               (1)  total energy of the orbit (alpha)

            \sstitem
               (2)  reference (osculating) epoch (t0)

            \sstitem
               (3-5)  position at reference epoch (r0)

            \sstitem
               (6-8)  velocity at reference epoch (v0)

            \sstitem
               (9)  heliocentric distance at reference epoch

            \sstitem
               (10)  r0.v0

            \sstitem
               (11)  date (t)

            \sstitem
               (12)  universal eccentric anomaly (psi) of date, approx
         }
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         Topocentric apparent RA (radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Topocentric apparent Dec (radians)
      }
      \sstsubsection{
         r = double $*$ (Returned)
      }{
         Distance from observer (AU)
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status: 0 = OK
         \sstitemlist{

            \sstitem
                  -1 = radius vector zero

            \sstitem
                  -2 = failed to converge
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          DATE is the instant for which the prediction is required.  It is
           in the TT timescale (formerly Ephemeris Time, ET) and is a
           Modified Julian Date (JD-2400000.5).

         \sstitem
          The longitude and latitude allow correction for geocentric
           parallax.  This is usually a small effect, but can become
           important for near-Earth asteroids.  Geocentric positions can be
           generated by appropriate use of routines palEpv (or palEvp) and
           palUe2pv.

         \sstitem
          The \texttt{"} universal\texttt{"}  elements are those which define the orbit for the
           purposes of the method of universal variables (see reference 2).
           They consist of the combined mass of the two bodies, an epoch,
           and the position and velocity vectors (arbitrary reference frame)
           at that epoch.  The parameter set used here includes also various
           quantities that can, in fact, be derived from the other
           information.  This approach is taken to avoiding unnecessary
           computation and loss of accuracy.  The supplementary quantities
           are (i) alpha, which is proportional to the total energy of the
           orbit, (ii) the heliocentric distance at epoch, (iii) the
           outwards component of the velocity at the given epoch, (iv) an
           estimate of psi, the \texttt{"} universal eccentric anomaly\texttt{"}  at a given
           date and (v) that date.

         \sstitem
          The universal elements are with respect to the J2000 equator and
           equinox.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          Sterne, Theodore E., \texttt{"} An Introduction to Celestial Mechanics\texttt{"} ,
           Interscience Publishers Inc., 1960.  Section 6.7, p199.

         \sstitem
          Everhart, E. \& Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
      }
   }
}
\sstroutine{
   palPm
}{
   Apply corrections for proper motion a star RA,Dec
}{
   \sstdescription{
      Apply corrections for proper motion to a star RA,Dec using the
      SOFA/ERFA routine eraStarpm.
   }
   \sstinvocation{
      void palPm ( double r0, double d0, double pr, double pd,
                   double px, double rv, double ep0, double ep1,
                   double $*$r1, double $*$d1 );
   }
   \sstarguments{
      \sstsubsection{
         r0 = double (Given)
      }{
         RA at epoch ep0 (radians)
      }
      \sstsubsection{
         d0 = double (Given)
      }{
         Dec at epoch ep0 (radians)
      }
      \sstsubsection{
         pr = double (Given)
      }{
         RA proper motion in radians per year.
      }
      \sstsubsection{
         pd = double (Given)
      }{
         Dec proper motion in radians per year.
      }
      \sstsubsection{
         px = double (Given)
      }{
         Parallax (arcsec)
      }
      \sstsubsection{
         rv = double (Given)
      }{
         Radial velocity (km/sec $+$ve if receding)
      }
      \sstsubsection{
         ep0 = double (Given)
      }{
         Start epoch in years, assumed to be Julian.
      }
      \sstsubsection{
         ep1 = double (Given)
      }{
         End epoch in years, assumed to be Julian.
      }
      \sstsubsection{
         r1 = double $*$ (Returned)
      }{
         RA at epoch ep1 (radians)
      }
      \sstsubsection{
         d1 = double $*$ (Returned)
      }{
         Dec at epoch ep1 (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraStarpm but ignores the status returns from that routine.
           In particular note that parallax should not be zero when the
           proper motions are non-zero. SLA/F allows parallax to be zero.

         \sstitem
          Assumes all epochs are Julian epochs.
      }
   }
}
\sstroutine{
   palPolmo
}{
   Correct for polar motion
}{
   \sstdescription{
      Polar motion:  correct site longitude and latitude for polar
      motion and calculate azimuth difference between celestial and
      terrestrial poles.
   }
   \sstinvocation{
      palPolmo ( double elongm, double phim, double xp, double yp,
                 double $*$elong, double $*$phi, double $*$daz );
   }
   \sstarguments{
      \sstsubsection{
         elongm = double (Given)
      }{
         Mean logitude of the observer (radians, east $+$ve)
      }
      \sstsubsection{
         phim = double (Given)
      }{
         Mean geodetic latitude of the observer (radians)
      }
      \sstsubsection{
         xp = double (Given)
      }{
         Polar motion x-coordinate (radians)
      }
      \sstsubsection{
         yp = double (Given)
      }{
         Polar motion y-coordinate (radians)
      }
      \sstsubsection{
         elong = double $*$ (Returned)
      }{
         True longitude of the observer (radians, east $+$ve)
      }
      \sstsubsection{
         phi = double $*$ (Returned)
      }{
         True geodetic latitude of the observer (radians)
      }
      \sstsubsection{
         daz = double $*$ (Returned)
      }{
         Azimuth correction (terrestrial-celestial, radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          \texttt{"} Mean\texttt{"}  longitude and latitude are the (fixed) values for the
           site\texttt{'} s location with respect to the IERS terrestrial reference
           frame;  the latitude is geodetic.  TAKE CARE WITH THE LONGITUDE
           SIGN CONVENTION.  The longitudes used by the present routine
           are east-positive, in accordance with geographical convention
           (and right-handed).  In particular, note that the longitudes
           returned by the sla\_OBS routine are west-positive, following
           astronomical usage, and must be reversed in sign before use in
           the present routine.

         \sstitem
          XP and YP are the (changing) coordinates of the Celestial
           Ephemeris Pole with respect to the IERS Reference Pole.
           XP is positive along the meridian at longitude 0 degrees,
           and YP is positive along the meridian at longitude
           270 degrees (i.e. 90 degrees west).  Values for XP,YP can
           be obtained from IERS circulars and equivalent publications;
           the maximum amplitude observed so far is about 0.3 arcseconds.

         \sstitem
          \texttt{"} True\texttt{"}  longitude and latitude are the (moving) values for
           the site\texttt{'} s location with respect to the celestial ephemeris
           pole and the meridian which corresponds to the Greenwich
           apparent sidereal time.  The true longitude and latitude
           link the terrestrial coordinates with the standard celestial
           models (for precession, nutation, sidereal time etc).

         \sstitem
          The azimuths produced by sla\_AOP and sla\_AOPQK are with
           respect to due north as defined by the Celestial Ephemeris
           Pole, and can therefore be called \texttt{"} celestial azimuths\texttt{"} .
           However, a telescope fixed to the Earth measures azimuth
           essentially with respect to due north as defined by the
           IERS Reference Pole, and can therefore be called \texttt{"} terrestrial
           azimuth\texttt{"} .  Uncorrected, this would manifest itself as a
           changing \texttt{"} azimuth zero-point error\texttt{"} .  The value DAZ is the
           correction to be added to a celestial azimuth to produce
           a terrestrial azimuth.

         \sstitem
          The present routine is rigorous.  For most practical
           purposes, the following simplified formulae provide an
           adequate approximation:

      }
        elong = elongm$+$xp$*$cos(elongm)-yp$*$sin(elongm)
        phi   = phim$+$(xp$*$sin(elongm)$+$yp$*$cos(elongm))$*$tan(phim)
        daz   = -sqrt(xp$*$xp$+$yp$*$yp)$*$cos(elongm-atan2(xp,yp))/cos(phim)

        An alternative formulation for DAZ is:

        x = cos(elongm)$*$cos(phim)
        y = sin(elongm)$*$cos(phim)
        daz = atan2(-x$*$yp-y$*$xp,x$*$x$+$y$*$y)

      \sstitemlist{

         \sstitem
          Reference:  Seidelmann, P.K. (ed), 1992.  \texttt{"} Explanatory Supplement
                       to the Astronomical Almanac\texttt{"} , ISBN 0-935702-68-7,
                       sections 3.27, 4.25, 4.52.
      }
   }
}
\sstroutine{
   palPrebn
}{
   Generate the matrix of precession between two objects (old)
}{
   \sstdescription{
      Generate the matrix of precession between two epochs,
      using the old, pre-IAU1976, Bessel-Newcomb model, using
      Kinoshita\texttt{'} s formulation
   }
   \sstinvocation{
      void palPrebn ( double bep0, double bep1, double rmatp[3][3] );
   }
   \sstarguments{
      \sstsubsection{
         bep0 = double (Given)
      }{
         Beginning Besselian epoch.
      }
      \sstsubsection{
         bep1 = double (Given)
      }{
         Ending Besselian epoch
      }
      \sstsubsection{
         rmatp = double[3][3] (Returned)
      }{
         precession matrix in the sense V(BEP1) = RMATP $*$ V(BEP0)
      }
   }
   \sstdiytopic{
      See Also
   }{
      Kinoshita, H. (1975) \texttt{'} Formulas for precession\texttt{'} , SAO Special
      Report No. 364, Smithsonian Institution Astrophysical
      Observatory, Cambridge, Massachusetts.
   }
}
\sstroutine{
   palPrec
}{
   Form the matrix of precession between two epochs (IAU 2006)
}{
   \sstdescription{
      The IAU 2006 precession matrix from ep0 to ep1 is found and
      returned. The matrix is in the sense  V(EP1)  =  RMATP $*$ V(EP0).
      The epochs are TDB (loosely TT) Julian epochs.

      Though the matrix method itself is rigorous, the precession
      angles are expressed through canonical polynomials which are
      valid only for a limited time span of a few hundred years around
      the current epoch.
   }
   \sstinvocation{
      palPrec( double ep0, double ep1, double rmatp[3][3] )
   }
   \sstarguments{
      \sstsubsection{
         ep0 = double (Given)
      }{
         Beginning epoch
      }
      \sstsubsection{
         ep1 = double (Given)
      }{
         Ending epoch
      }
      \sstsubsection{
         rmatp = double[3][3] (Returned)
      }{
         Precession matrix
      }
   }
}
\sstroutine{
   palPreces
}{
   Precession - either FK4 or FK5 as required
}{
   \sstdescription{
      Precess coordinates using the appropriate system and epochs.
   }
   \sstinvocation{
      void palPreces ( const char sys[3], double ep0, double ep1,
                       double $*$ra, double $*$dc );
   }
   \sstarguments{
      \sstsubsection{
         sys = const char [3] (Given)
      }{
         Precession to be applied: FK4 or FK5. Case insensitive.
      }
      \sstsubsection{
         ep0 = double (Given)
      }{
         Starting epoch.
      }
      \sstsubsection{
         ep1 = double (Given)
      }{
         Ending epoch
      }
      \sstsubsection{
         ra = double $*$ (Given \& Returned)
      }{
         On input the RA mean equator \& equinox at epoch ep0. On exit
         the RA mean equator \& equinox of epoch ep1.
      }
      \sstsubsection{
         dec = double $*$ (Given \& Returned)
      }{
         On input the dec mean equator \& equinox at epoch ep0. On exit
         the dec mean equator \& equinox of epoch ep1.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses palPrec for FK5 data and palPrebn for FK4 data.

         \sstitem
          The epochs are Besselian if SYSTEM=\texttt{'} FK4\texttt{'}  and Julian if \texttt{'} FK5\texttt{'} .
            For example, to precess coordinates in the old system from
            equinox 1900.0 to 1950.0 the call would be:
                 palPreces( \texttt{"} FK4\texttt{"} , 1900.0, 1950.0, \&ra, \&dc );

         \sstitem
          This routine will NOT correctly convert between the old and
           the new systems - for example conversion from B1950 to J2000.
           For these purposes see palFk425, palFk524, palFk45z and
           palFk54z.

         \sstitem
          If an invalid SYSTEM is supplied, values of -99D0,-99D0 will
           be returned for both RA and DC.
      }
   }
}
\sstroutine{
   palPrenut
}{
   Form the matrix of bias-precession-nutation (IAU 2006/2000A)
}{
   \sstdescription{
      Form the matrix of bias-precession-nutation (IAU 2006/2000A).
      The epoch and date are TT (but TDB is usually close enough).
      The matrix is in the sense   v(true)  =  rmatpn $*$ v(mean).
   }
   \sstinvocation{
      void palPrenut( double epoch, double date, double rmatpn[3][3] )
   }
   \sstarguments{
      \sstsubsection{
         epoch = double (Returned)
      }{
         Julian epoch for mean coordinates.
      }
      \sstsubsection{
         date = double (Returned)
      }{
         Modified Julian Date (JD-2400000.5) for true coordinates.
      }
      \sstsubsection{
         rmatpn = double[3][3] (Returned)
      }{
         combined NPB matrix
      }
   }
}
\sstroutine{
   palPv2el
}{
   Position velocity to heliocentirc osculating elements
}{
   \sstdescription{
      Heliocentric osculating elements obtained from instantaneous position
      and velocity.
   }
   \sstinvocation{
      void palPv2el ( const double pv[6], double date, double pmass, int jformr,
                      int $*$jform, double $*$epoch, double $*$orbinc,
                      double $*$anode, double $*$perih, double $*$aorq, double $*$e,
                      double $*$aorl, double $*$dm, int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         pv = const double [6] (Given)
      }{
         Heliocentric x,y,z,xdot,ydot,zdot of date,
         J2000 equatorial triad (AU,AU/s; Note 1)
      }
      \sstsubsection{
         date = double (Given)
      }{
         Date (TT Modified Julian Date = JD-2400000.5)
      }
      \sstsubsection{
         pmass = double (Given)
      }{
         Mass of the planet (Sun=1; Note 2)
      }
      \sstsubsection{
         jformr = int (Given)
      }{
         Requested element set (1-3; Note 3)
      }
      \sstsubsection{
         jform = int $*$ (Returned)
      }{
         Element set actually returned (1-3; Note 4)
      }
      \sstsubsection{
         epoch = double $*$ (Returned)
      }{
         Epoch of elements (TT MJD)
      }
      \sstsubsection{
         orbinc = double $*$ (Returned)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode = double $*$ (Returned)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih = double $*$ (Returned)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq = double $*$ (Returned)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e = double $*$ (Returned)
      }{
         eccentricity
      }
      \sstsubsection{
         aorl = double $*$ (Returned)
      }{
         mean anomaly or longitude (radians, JFORM=1,2 only)
      }
      \sstsubsection{
         dm = double $*$ (Returned)
      }{
         daily motion (radians, JFORM=1 only)
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:  0 = OK
         \sstitemlist{

            \sstitem
                    -1 = illegal PMASS

            \sstitem
                    -2 = illegal JFORMR

            \sstitem
                    -3 = position/velocity out of range
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The PV 6-vector is with respect to the mean equator and equinox of
           epoch J2000.  The orbital elements produced are with respect to
           the J2000 ecliptic and mean equinox.

         \sstitem
          The mass, PMASS, is important only for the larger planets.  For
           most purposes (e.g. asteroids) use 0D0.  Values less than zero
           are illegal.

         \sstitem
          Three different element-format options are supported:

      }
        Option JFORM=1, suitable for the major planets:

        EPOCH  = epoch of elements (TT MJD)
        ORBINC = inclination i (radians)
        ANODE  = longitude of the ascending node, big omega (radians)
        PERIH  = longitude of perihelion, curly pi (radians)
        AORQ   = mean distance, a (AU)
        E      = eccentricity, e
        AORL   = mean longitude L (radians)
        DM     = daily motion (radians)

        Option JFORM=2, suitable for minor planets:

        EPOCH  = epoch of elements (TT MJD)
        ORBINC = inclination i (radians)
        ANODE  = longitude of the ascending node, big omega (radians)
        PERIH  = argument of perihelion, little omega (radians)
        AORQ   = mean distance, a (AU)
        E      = eccentricity, e
        AORL   = mean anomaly M (radians)

        Option JFORM=3, suitable for comets:

        EPOCH  = epoch of perihelion (TT MJD)
        ORBINC = inclination i (radians)
        ANODE  = longitude of the ascending node, big omega (radians)
        PERIH  = argument of perihelion, little omega (radians)
        AORQ   = perihelion distance, q (AU)
        E      = eccentricity, e

      \sstitemlist{

         \sstitem
          It may not be possible to generate elements in the form
           requested through JFORMR.  The caller is notified of the form
           of elements actually returned by means of the JFORM argument:

      }
         JFORMR   JFORM     meaning

           1        1       OK - elements are in the requested format
           1        2       never happens
           1        3       orbit not elliptical

           2        1       never happens
           2        2       OK - elements are in the requested format
           2        3       orbit not elliptical

           3        1       never happens
           3        2       never happens
           3        3       OK - elements are in the requested format

      \sstitemlist{

         \sstitem
          The arguments returned for each value of JFORM (cf Note 5: JFORM
           may not be the same as JFORMR) are as follows:

      }
          JFORM         1              2              3
          EPOCH         t0             t0             T
          ORBINC        i              i              i
          ANODE         Omega          Omega          Omega
          PERIH         curly pi       omega          omega
          AORQ          a              a              q
          E             e              e              e
          AORL          L              M              -
          DM            n              -              -

        where:

          t0           is the epoch of the elements (MJD, TT)
          T              \texttt{"}     epoch of perihelion (MJD, TT)
          i              \texttt{"}     inclination (radians)
          Omega          \texttt{"}     longitude of the ascending node (radians)
          curly pi       \texttt{"}     longitude of perihelion (radians)
          omega          \texttt{"}     argument of perihelion (radians)
          a              \texttt{"}     mean distance (AU)
          q              \texttt{"}     perihelion distance (AU)
          e              \texttt{"}     eccentricity
          L              \texttt{"}     longitude (radians, 0-2pi)
          M              \texttt{"}     mean anomaly (radians, 0-2pi)
          n              \texttt{"}     daily motion (radians)
      \sstitemlist{

         \sstitem
              means no value is set

         \sstitem
          At very small inclinations, the longitude of the ascending node
           ANODE becomes indeterminate and under some circumstances may be
           set arbitrarily to zero.  Similarly, if the orbit is close to
           circular, the true anomaly becomes indeterminate and under some
           circumstances may be set arbitrarily to zero.  In such cases,
           the other elements are automatically adjusted to compensate,
           and so the elements remain a valid description of the orbit.

         \sstitem
          The osculating epoch for the returned elements is the argument
           DATE.

         \sstitem
          Reference:  Sterne, Theodore E., \texttt{"} An Introduction to Celestial
                       Mechanics\texttt{"} , Interscience Publishers, 1960
      }
   }
}
\sstroutine{
   palPv2ue
}{
   Universal elements to position and velocity
}{
   \sstdescription{
      Construct a universal element set based on an instantaneous position
      and velocity.
   }
   \sstinvocation{
      void palPv2ue( const double pv[6], double date, double pmass,
                     double u[13], int $*$ jstat );
   }
   \sstarguments{
      \sstsubsection{
         pv = double [6] (Given)
      }{
         Heliocentric x,y,z,xdot,ydot,zdot of date, (AU,AU/s; Note 1)
      }
      \sstsubsection{
         date = double (Given)
      }{
         Date (TT modified Julian Date = JD-2400000.5)
      }
      \sstsubsection{
         pmass = double (Given)
      }{
         Mass of the planet (Sun=1; note 2)
      }
      \sstsubsection{
         u = double [13] (Returned)
      }{
         Universal orbital elements (Note 3)

         \sstitemlist{

            \sstitem
               (0)  combined mass (M$+$m)

            \sstitem
               (1)  total energy of the orbit (alpha)

            \sstitem
               (2)  reference (osculating) epoch (t0)

            \sstitem
               (3-5)  position at reference epoch (r0)

            \sstitem
               (6-8)  velocity at reference epoch (v0)

            \sstitem
               (9)  heliocentric distance at reference epoch

            \sstitem
               (10)  r0.v0

            \sstitem
               (11)  date (t)

            \sstitem
               (12)  universal eccentric anomaly (psi) of date, approx
         }
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status: 0 = OK
         \sstitemlist{

            \sstitem
                    -1 = illegal PMASS

            \sstitem
                    -2 = too close to Sun

            \sstitem
                    -3 = too slow
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The PV 6-vector can be with respect to any chosen inertial frame,
           and the resulting universal-element set will be with respect to
           the same frame.  A common choice will be mean equator and ecliptic
           of epoch J2000.

         \sstitem
          The mass, PMASS, is important only for the larger planets.  For
           most purposes (e.g. asteroids) use 0D0.  Values less than zero
           are illegal.

         \sstitem
          The \texttt{"} universal\texttt{"}  elements are those which define the orbit for the
           purposes of the method of universal variables (see reference).
           They consist of the combined mass of the two bodies, an epoch,
           and the position and velocity vectors (arbitrary reference frame)
           at that epoch.  The parameter set used here includes also various
           quantities that can, in fact, be derived from the other
           information.  This approach is taken to avoiding unnecessary
           computation and loss of accuracy.  The supplementary quantities
           are (i) alpha, which is proportional to the total energy of the
           orbit, (ii) the heliocentric distance at epoch, (iii) the
           outwards component of the velocity at the given epoch, (iv) an
           estimate of psi, the \texttt{"} universal eccentric anomaly\texttt{"}  at a given
           date and (v) that date.

         \sstitem
          Reference:  Everhart, E. \& Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
      }
   }
}
\sstroutine{
   palPvobs
}{
   Position and velocity of an observing station
}{
   \sstdescription{
      Returns the position and velocity of an observing station.
   }
   \sstinvocation{
      palPvobs( double p, double h, double stl, double pv[6] )
   }
   \sstarguments{
      \sstsubsection{
         p = double (Given)
      }{
         Latitude (geodetic, radians).
      }
      \sstsubsection{
         h = double (Given)
      }{
         Height above reference spheroid (geodetic, metres).
      }
      \sstsubsection{
         stl = double (Given)
      }{
         Local apparent sidereal time (radians).
      }
      \sstsubsection{
         pv = double[ 6 ] (Returned)
      }{
         position/velocity 6-vector (AU, AU/s, true equator
                                     and equinox of date).
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The WGS84 reference ellipsoid is used.
      }
   }
}
\sstroutine{
   palPxy
}{
   Given arrays of \texttt{"} expected\texttt{"}  and \texttt{"} measured\texttt{"}  [x,y] coordinates, and a
   linear model relating them (as produced by palFitxy), compute
   the array of \texttt{"} predicted\texttt{"}  coordinates and the RMS residuals
}{
   \sstdescription{
      The model is supplied in the array coeffs.  Naming the
      elements of coeffs as follows:
      \newline
      \newline
      \hspace*{1.5 em}
         coeffs[0] = A
      \newline
      \hspace*{1.5 em}
         coeffs[1] = B
      \newline
      \hspace*{1.5 em}
         coeffs[2] = C
      \newline
      \hspace*{1.5 em}
         coeffs[3] = D
      \newline
      \hspace*{1.5 em}
         coeffs[4] = E
      \newline
      \hspace*{1.5 em}
         coeffs[5] = F
      \newline
      \newline
      the model is applied thus:
      \newline
      \newline
      \hspace*{1.5 em}
         xp = A $+$ B $*$ xm $+$ C $*$ ym
      \newline
      \hspace*{1.5 em}
         yp = D $+$ E $*$ xm $+$ F $*$ ym
      \newline
      \newline
      The residuals are (xp - xe) and (yp - ye).

      If np is less than or equal to zero, no coordinates are
      transformed, and the RMS residuals are all zero.
   }
   \sstinvocation{
      palPxy ( int np, double xye[][2], double xym[][2], double coeffs[6],
               double xyp[][2], double $*$xrms, double $*$yrms, double $*$rrms )
   }
   \sstarguments{
      \sstsubsection{
         np = int (Given)
      }{
         number of samples
      }
      \sstsubsection{
         xye = double[np][2] (Given)
      }{
         expected [x,y] for each sample
      }
      \sstsubsection{
         xym = double[np][2] (Given)
      }{
         measured [x,y] for each sample
      }
      \sstsubsection{
         coeffs = double[6]
      }{
         coefficients of model (see below)
      }
      \sstsubsection{
         xyp = double[np][2] (Returned)
      }{
         predicted [x,y] for each sample
      }
      \sstsubsection{
         xrms = double (Returned)
      }{
         RMS in x
      }
      \sstsubsection{
         yrms = double (Returned)
      }{
         RMS in y
      }
      \sstsubsection{
         rrms = double (Returned)
      }{
         total RMS (vector sum of xrms and yrms)
      }
   }
   \sstdiytopic{
      See also
   }{
      palFitxy, palInvf, palXy2xy and palDcmpf
   }
}
\sstroutine{
   palRanorm
}{
   Reduce angle to be within [0, 2$*$pi) (single precision)
}{
   \sstdescription{
      The result is \texttt{"} angle\texttt{"}  expressed in the range 0 to 2$*$pi.
   }
   \sstinvocation{
      float palRanorm ( float angle );
   }
   \sstarguments{
      \sstsubsection{
         angle = float (Given)
      }{
         Angle to be reduced (radians)
      }
   }
}
\sstroutine{
   palRdplan
}{
   Approximate topocentric apparent RA,Dec of a planet
}{
   \sstdescription{
      Approximate topocentric apparent RA,Dec of a planet, and its
      angular diameter.
   }
   \sstinvocation{
      void palRdplan( double date, int np, double elong, double phi,
                      double $*$ ra, double $*$ dec, double $*$ diam );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         MJD of observation (JD-2400000.5) in TDB. For all practical
         purposes TT can be used instead of TDB, and for many applications
         UT will do (except for the Moon).
      }
      \sstsubsection{
         np = int (Given)
      }{
         Planet: 1 = Mercury
                 2 = Venus
                 3 = Moon
                 4 = Mars
                 5 = Jupiter
                 6 = Saturn
                 7 = Uranus
                 8 = Neptune
              else = Sun
      }
      \sstsubsection{
         elong = double (Given)
      }{
         Observer\texttt{'} s east longitude (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observer\texttt{'} s geodetic latitude (radians)
      }
      \sstsubsection{
         ra = double $*$ (Returned)
      }{
         RA (topocentric apparent, radians)
      }
      \sstsubsection{
         dec = double $*$ (Returned)
      }{
         Dec (topocentric apparent, radians)
      }
      \sstsubsection{
         diam = double $*$ (Returned)
      }{
         Angular diameter (equatorial, radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Unlike with slaRdplan, Pluto is not supported.

         \sstitem
          The longitude and latitude allow correction for geocentric
           parallax.  This is a major effect for the Moon, but in the
           context of the limited accuracy of the present routine its
           effect on planetary positions is small (negligible for the
           outer planets).  Geocentric positions can be generated by
           appropriate use of the routines palDmoon and eraPlan94.
      }
   }
}
\sstroutine{
   palRefco
}{
   Determine constants in atmospheric refraction model
}{
   \sstdescription{
      Determine the constants A and B in the atmospheric refraction
      model dZ = A tan Z $+$ B tan$*$$*$3 Z.

      Z is the \texttt{"} observed\texttt{"}  zenith distance (i.e. affected by refraction)
      and dZ is what to add to Z to give the \texttt{"} topocentric\texttt{"}  (i.e. in vacuo)
      zenith distance.
   }
   \sstinvocation{
      void palRefco ( double hm, double tdk, double pmb, double rh,
                      double wl, double phi, double tlr, double eps,
                      double $*$refa, double $*$refb );
   }
   \sstarguments{
      \sstsubsection{
         hm = double (Given)
      }{
         Height of the observer above sea level (metre)
      }
      \sstsubsection{
         tdk = double (Given)
      }{
         Ambient temperature at the observer (K)
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         Pressure at the observer (millibar)
      }
      \sstsubsection{
         rh = double (Given)
      }{
         Relative humidity at the observer (range 0-1)
      }
      \sstsubsection{
         wl = double (Given)
      }{
         Effective wavelength of the source (micrometre)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Latitude of the observer (radian, astronomical)
      }
      \sstsubsection{
         tlr = double (Given)
      }{
         Temperature lapse rate in the troposphere (K/metre)
      }
      \sstsubsection{
         eps = double (Given)
      }{
         Precision required to terminate iteration (radian)
      }
      \sstsubsection{
         refa = double $*$ (Returned)
      }{
         tan Z coefficient (radian)
      }
      \sstsubsection{
         refb = double $*$ (Returned)
      }{
         tan$*$$*$3 Z coefficient (radian)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Typical values for the TLR and EPS arguments might be 0.0065 and
         1E-10 respectively.

         \sstitem
          The radio refraction is chosen by specifying WL $>$ 100 micrometres.

         \sstitem
          The routine is a slower but more accurate alternative to the
         palRefcoq routine.  The constants it produces give perfect
         agreement with palRefro at zenith distances arctan(1) (45 deg)
         and arctan(4) (about 76 deg).  It achieves 0.5 arcsec accuracy
         for ZD $<$ 80 deg, 0.01 arcsec accuracy for ZD $<$ 60 deg, and
         0.001 arcsec accuracy for ZD $<$ 45 deg.
      }
   }
}
\sstroutine{
   palRefro
}{
   Atmospheric refraction for radio and optical/IR wavelengths
}{
   \sstdescription{
      Calculates the atmospheric refraction for radio and optical/IR
      wavelengths.
   }
   \sstinvocation{
      void palRefro( double zobs, double hm, double tdk, double pmb,
                     double rh, double wl, double phi, double tlr,
                     double eps, double $*$ ref ) \{
   }
   \sstarguments{
      \sstsubsection{
         zobs = double (Given)
      }{
         Observed zenith distance of the source (radian)
      }
      \sstsubsection{
         hm = double (Given)
      }{
         Height of the observer above sea level (metre)
      }
      \sstsubsection{
         tdk = double (Given)
      }{
         Ambient temperature at the observer (K)
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         Pressure at the observer (millibar)
      }
      \sstsubsection{
         rh = double (Given)
      }{
         Relative humidity at the observer (range 0-1)
      }
      \sstsubsection{
         wl = double (Given)
      }{
         Effective wavelength of the source (micrometre)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Latitude of the observer (radian, astronomical)
      }
      \sstsubsection{
         tlr = double (Given)
      }{
         Temperature lapse rate in the troposphere (K/metre)
      }
      \sstsubsection{
         eps = double (Given)
      }{
         Precision required to terminate iteration (radian)
      }
      \sstsubsection{
         ref = double $*$ (Returned)
      }{
         Refraction: in vacuao ZD minus observed ZD (radian)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          A suggested value for the TLR argument is 0.0065.  The
         refraction is significantly affected by TLR, and if studies
         of the local atmosphere have been carried out a better TLR
         value may be available.  The sign of the supplied TLR value
         is ignored.

         \sstitem
          A suggested value for the EPS argument is 1E-8.  The result is
         usually at least two orders of magnitude more computationally
         precise than the supplied EPS value.

         \sstitem
          The routine computes the refraction for zenith distances up
         to and a little beyond 90 deg using the method of Hohenkerk
         and Sinclair (NAO Technical Notes 59 and 63, subsequently adopted
         in the Explanatory Supplement, 1992 edition - see section 3.281).

         \sstitem
          The code is a development of the optical/IR refraction subroutine
         AREF of C.Hohenkerk (HMNAO, September 1984), with extensions to
         support the radio case.  Apart from merely cosmetic changes, the
         following modifications to the original HMNAO optical/IR refraction
         code have been made:

      }
      .  The angle arguments have been changed to radians.

      .  Any value of ZOBS is allowed (see note 6, below).

      .  Other argument values have been limited to safe values.

      .  Murray\texttt{'} s values for the gas constants have been used
         (Vectorial Astrometry, Adam Hilger, 1983).

      .  The numerical integration phase has been rearranged for
         extra clarity.

      .  A better model for Ps(T) has been adopted (taken from
         Gill, Atmosphere-Ocean Dynamics, Academic Press, 1982).

      .  More accurate expressions for Pwo have been adopted
         (again from Gill 1982).

      .  The formula for the water vapour pressure, given the
         saturation pressure and the relative humidity, is from
         Crane (1976), expression 2.5.5.

      .  Provision for radio wavelengths has been added using
         expressions devised by A.T.Sinclair, RGO (private
         communication 1989).  The refractivity model currently
         used is from J.M.Rueger, \texttt{"} Refractive Index Formulae for
         Electronic Distance Measurement with Radio and Millimetre
         Waves\texttt{"} , in Unisurv Report S-68 (2002), School of Surveying
         and Spatial Information Systems, University of New South
         Wales, Sydney, Australia.

      .  The optical refractivity for dry air is from Resolution 3 of
         the International Association of Geodesy adopted at the XXIIth
         General Assembly in Birmingham, UK, 1999.

      .  Various small changes have been made to gain speed.

      \sstitemlist{

         \sstitem
          The radio refraction is chosen by specifying WL $>$ 100 micrometres.
         Because the algorithm takes no account of the ionosphere, the
         accuracy deteriorates at low frequencies, below about 30 MHz.

         \sstitem
          Before use, the value of ZOBS is expressed in the range $+$/- pi.
         If this ranged ZOBS is -ve, the result REF is computed from its
         absolute value before being made -ve to match.  In addition, if
         it has an absolute value greater than 93 deg, a fixed REF value
         equal to the result for ZOBS = 93 deg is returned, appropriately
         signed.

         \sstitem
          As in the original Hohenkerk and Sinclair algorithm, fixed values
         of the water vapour polytrope exponent, the height of the
         tropopause, and the height at which refraction is negligible are
         used.

         \sstitem
          The radio refraction has been tested against work done by
         Iain Coulson, JACH, (private communication 1995) for the
         James Clerk Maxwell Telescope, Mauna Kea.  For typical conditions,
         agreement at the 0.1 arcsec level is achieved for moderate ZD,
         worsening to perhaps 0.5-1.0 arcsec at ZD 80 deg.  At hot and
         humid sea-level sites the accuracy will not be as good.

         \sstitem
          It should be noted that the relative humidity RH is formally
         defined in terms of \texttt{"} mixing ratio\texttt{"}  rather than pressures or
         densities as is often stated.  It is the mass of water per unit
         mass of dry air divided by that for saturated air at the same
         temperature and pressure (see Gill 1982).

         \sstitem
          The algorithm is designed for observers in the troposphere. The
         supplied temperature, pressure and lapse rate are assumed to be
         for a point in the troposphere and are used to define a model
         atmosphere with the tropopause at 11km altitude and a constant
         temperature above that.  However, in practice, the refraction
         values returned for stratospheric observers, at altitudes up to
         25km, are quite usable.
      }
   }
}
\sstroutine{
   palRefv
}{
   Adjust an unrefracted Cartesian vector to include the effect of atmospheric refraction
}{
   \sstdescription{
      Adjust an unrefracted Cartesian vector to include the effect of
      atmospheric refraction, using the simple A tan Z $+$ B tan$*$$*$3 Z
      model.
   }
   \sstinvocation{
      void palRefv ( double vu[3], double refa, double refb, double vr[3] );
   }
   \sstarguments{
      \sstsubsection{
         vu[3] = double (Given)
      }{
         Unrefracted position of the source (Az/El 3-vector)
      }
      \sstsubsection{
         refa = double (Given)
      }{
         tan Z coefficient (radian)
      }
      \sstsubsection{
         refb = double (Given)
      }{
         tan$*$$*$3 Z coefficient (radian)
      }
      \sstsubsection{
         vr[3] = double (Returned)
      }{
         Refracted position of the source (Az/El 3-vector)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          This routine applies the adjustment for refraction in the
         opposite sense to the usual one - it takes an unrefracted
         (in vacuo) position and produces an observed (refracted)
         position, whereas the A tan Z $+$ B tan$*$$*$3 Z model strictly
         applies to the case where an observed position is to have the
         refraction removed.  The unrefracted to refracted case is
         harder, and requires an inverted form of the text-book
         refraction models;  the algorithm used here is equivalent to
         one iteration of the Newton-Raphson method applied to the above
         formula.

         \sstitem
          Though optimized for speed rather than precision, the present
         routine achieves consistency with the refracted-to-unrefracted
         A tan Z $+$ B tan$*$$*$3 Z model at better than 1 microarcsecond within
         30 degrees of the zenith and remains within 1 milliarcsecond to
         beyond ZD 70 degrees.  The inherent accuracy of the model is, of
         course, far worse than this - see the documentation for palRefco
         for more information.

         \sstitem
          At low elevations (below about 3 degrees) the refraction
         correction is held back to prevent arithmetic problems and
         wildly wrong results.  For optical/IR wavelengths, over a wide
         range of observer heights and corresponding temperatures and
         pressures, the following levels of accuracy (arcsec, worst case)
         are achieved, relative to numerical integration through a model
         atmosphere:

      }
               ZD    error

               80      0.7
               81      1.3
               82      2.5
               83      5
               84     10
               85     20
               86     55
               87    160
               88    360
               89    640
               90   1100
               91   1700         \} relevant only to
               92   2600         \} high-elevation sites

      The results for radio are slightly worse over most of the range,
      becoming significantly worse below ZD=88 and unusable beyond
      ZD=90.

      \sstitemlist{

         \sstitem
          See also the routine palRefz, which performs the adjustment to
         the zenith distance rather than in Cartesian Az/El coordinates.
         The present routine is faster than palRefz and, except very low down,
         is equally accurate for all practical purposes.  However, beyond
         about ZD 84 degrees palRefz should be used, and for the utmost
         accuracy iterative use of palRefro should be considered.
      }
   }
}
\sstroutine{
   palRefz
}{
   Adjust unrefracted zenith distance
}{
   \sstdescription{
      Adjust an unrefracted zenith distance to include the effect of
      atmospheric refraction, using the simple A tan Z $+$ B tan$*$$*$3 Z
      model (plus special handling for large ZDs).
   }
   \sstinvocation{
      void palRefz ( double zu, double refa, double refb, double $*$zr );
   }
   \sstarguments{
      \sstsubsection{
         zu = double (Given)
      }{
         Unrefracted zenith distance of the source (radians)
      }
      \sstsubsection{
         refa = double (Given)
      }{
         tan Z coefficient (radians)
      }
      \sstsubsection{
         refb = double (Given)
      }{
         tan$*$$*$3 Z coefficient (radian)
      }
      \sstsubsection{
         zr = double $*$ (Returned)
      }{
         Refracted zenith distance (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          This routine applies the adjustment for refraction in the
         opposite sense to the usual one - it takes an unrefracted
         (in vacuo) position and produces an observed (refracted)
         position, whereas the A tan Z $+$ B tan$*$$*$3 Z model strictly
         applies to the case where an observed position is to have the
         refraction removed.  The unrefracted to refracted case is
         harder, and requires an inverted form of the text-book
         refraction models;  the formula used here is based on the
         Newton-Raphson method.  For the utmost numerical consistency
         with the refracted to unrefracted model, two iterations are
         carried out, achieving agreement at the 1D-11 arcseconds level
         for a ZD of 80 degrees.  The inherent accuracy of the model
         is, of course, far worse than this - see the documentation for
         palRefco for more information.

         \sstitem
          At ZD 83 degrees, the rapidly-worsening A tan Z $+$ B tan$\wedge$3 Z
         model is abandoned and an empirical formula takes over.  For
         optical/IR wavelengths, over a wide range of observer heights and
         corresponding temperatures and pressures, the following levels of
         accuracy (arcsec, worst case) are achieved, relative to numerical
         integration through a model atmosphere:

      }
               ZR    error

               80      0.7
               81      1.3
               82      2.4
               83      4.7
               84      6.2
               85      6.4
               86      8
               87     10
               88     15
               89     30
               90     60
               91    150         \} relevant only to
               92    400         \} high-elevation sites

      For radio wavelengths the errors are typically 50\% larger than
      the optical figures and by ZD 85 deg are twice as bad, worsening
      rapidly below that.  To maintain 1 arcsec accuracy down to ZD=85
      at the Green Bank site, Condon (2004) has suggested amplifying
      the amount of refraction predicted by palRefz below 10.8 deg
      elevation by the factor (1$+$0.00195$*$(10.8-E\_t)), where E\_t is the
      unrefracted elevation in degrees.

      The high-ZD model is scaled to match the normal model at the
      transition point;  there is no glitch.

      \sstitemlist{

         \sstitem
          Beyond 93 deg zenith distance, the refraction is held at its
         93 deg value.

         \sstitem
          See also the routine palRefv, which performs the adjustment in
         Cartesian Az/El coordinates, and with the emphasis on speed
         rather than numerical accuracy.
      }
   }
   \sstdiytopic{
      References
   }{
      Condon,J.J., Refraction Corrections for the GBT, PTCS/PN/35.2,
      NRAO Green Bank, 2004.
   }
}
\sstroutine{
   palRverot
}{
   Velocity component in a given direction due to Earth rotation
}{
   \sstdescription{
      Calculate the velocity component in a given direction due to Earth
      rotation.

      The simple algorithm used assumes a spherical Earth, of
      a radius chosen to give results accurate to about 0.0005 km/s
      for observing stations at typical latitudes and heights.  For
      applications requiring greater precision, use the routine
      palPvobs.
   }
   \sstinvocation{
      double palRverot ( double phi, double ra, double da, double st );
   }
   \sstarguments{
      \sstsubsection{
         phi = double (Given)
      }{
         latitude of observing station (geodetic) (radians)
      }
      \sstsubsection{
         ra = double (Given)
      }{
         apparent RA (radians)
      }
      \sstsubsection{
         da = double (Given)
      }{
         apparent Dec (radians)
      }
      \sstsubsection{
         st = double (Given)
      }{
         Local apparent sidereal time.
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         palRverot = double
      }{
         Component of Earth rotation in direction RA,DA (km/s).
         The result is $+$ve when the observatory is receding from the
         given point on the sky.
      }
   }
}
\sstroutine{
   palRvgalc
}{
   Velocity component in a given direction due to the rotation
   of the Galaxy
}{
   \sstdescription{
      This function returns the Component of dynamical LSR motion in
      the direction of R2000,D2000. The result is $+$ve when the dynamical
      LSR is receding from the given point on the sky.
   }
   \sstinvocation{
      double palRvgalc( double r2000, double d2000 )
   }
   \sstarguments{
      \sstsubsection{
         r2000 = double (Given)
      }{
         J2000.0 mean RA (radians)
      }
      \sstsubsection{
         d2000 = double (Given)
      }{
         J2000.0 mean Dec (radians)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Component of dynamical LSR motion in direction R2000,D2000 (km/s).
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The Local Standard of Rest used here is a point in the
         vicinity of the Sun which is in a circular orbit around
         the Galactic centre.  Sometimes called the \texttt{"} dynamical\texttt{"}  LSR,
         it is not to be confused with a \texttt{"} kinematical\texttt{"}  LSR, which
         is the mean standard of rest of star catalogues or stellar
         populations.
      }
   }
   \sstdiytopic{
      Reference
   }{
      \sstitemlist{

         \sstitem
          The orbital speed of 220 km/s used here comes from Kerr \&
         Lynden-Bell (1986), MNRAS, 221, p1023.
      }
   }
}
\sstroutine{
   palRvlg
}{
   Velocity component in a given direction due to Galactic rotation
   and motion of the local group
}{
   \sstdescription{
      This function returns the velocity component in a given
      direction due to the combination of the rotation of the
      Galaxy and the motion of the Galaxy relative to the mean
      motion of the local group. The result is $+$ve when the Sun
      is receding from the given point on the sky.
   }
   \sstinvocation{
      double palRvlg( double r2000, double d2000 )
   }
   \sstarguments{
      \sstsubsection{
         r2000 = double (Given)
      }{
         J2000.0 mean RA (radians)
      }
      \sstsubsection{
         d2000 = double (Given)
      }{
         J2000.0 mean Dec (radians)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Component of SOLAR motion in direction R2000,D2000 (km/s).
      }{
      }
   }
   \sstdiytopic{
      Reference
   }{
      \sstitemlist{

         \sstitem
          IAU Trans 1976, 168, p201.
      }
   }
}
\sstroutine{
   palRvlsrd
}{
   Velocity component in a given direction due to the Sun\texttt{'} s motion
   with respect to the dynamical Local Standard of Rest
}{
   \sstdescription{
      This function returns the velocity component in a given direction
      due to the Sun\texttt{'} s motion with respect to the dynamical Local Standard
      of Rest. The result is $+$ve when the Sun is receding from the given
      point on the sky.
   }
   \sstinvocation{
      double palRvlsrd( double r2000, double d2000 )
   }
   \sstarguments{
      \sstsubsection{
         r2000 = double (Given)
      }{
         J2000.0 mean RA (radians)
      }
      \sstsubsection{
         d2000 = double (Given)
      }{
         J2000.0 mean Dec (radians)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Component of \texttt{"} peculiar\texttt{"}  solar motion in direction R2000,D2000 (km/s).
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The Local Standard of Rest used here is the \texttt{"} dynamical\texttt{"}  LSR,
         a point in the vicinity of the Sun which is in a circular orbit
         around the Galactic centre.  The Sun\texttt{'} s motion with respect to the
         dynamical LSR is called the \texttt{"} peculiar\texttt{"}  solar motion.

         \sstitem
          There is another type of LSR, called a \texttt{"} kinematical\texttt{"}  LSR.  A
         kinematical LSR is the mean standard of rest of specified star
         catalogues or stellar populations, and several slightly different
         kinematical LSRs are in use.  The Sun\texttt{'} s motion with respect to an
         agreed kinematical LSR is known as the \texttt{"} standard\texttt{"}  solar motion.
         To obtain a radial velocity correction with respect to an adopted
         kinematical LSR use the routine palRvlsrk.
      }
   }
   \sstdiytopic{
      Reference
   }{
      \sstitemlist{

         \sstitem
          Delhaye (1965), in \texttt{"} Stars and Stellar Systems\texttt{"} , vol 5, p73.
      }
   }
}
\sstroutine{
   palRvlsrk
}{
   Velocity component in a given direction due to the Sun\texttt{'} s motion
   with respect to an adopted kinematic Local Standard of Rest
}{
   \sstdescription{
      This function returns the velocity component in a given direction
      due to the Sun\texttt{'} s motion with respect to an adopted kinematic
      Local Standard of Rest. The result is $+$ve when the Sun is receding
      from the given point on the sky.
   }
   \sstinvocation{
      double palRvlsrk( double r2000, double d2000 )
   }
   \sstarguments{
      \sstsubsection{
         r2000 = double (Given)
      }{
         J2000.0 mean RA (radians)
      }
      \sstsubsection{
         d2000 = double (Given)
      }{
         J2000.0 mean Dec (radians)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Component of \texttt{"} standard\texttt{"}  solar motion in direction R2000,D2000 (km/s).
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The Local Standard of Rest used here is one of several
         \texttt{"} kinematical\texttt{"}  LSRs in common use.  A kinematical LSR is the mean
         standard of rest of specified star catalogues or stellar
         populations.  The Sun\texttt{'} s motion with respect to a kinematical LSR
         is known as the \texttt{"} standard\texttt{"}  solar motion.

         \sstitem
          There is another sort of LSR, the \texttt{"} dynamical\texttt{"}  LSR, which is a
         point in the vicinity of the Sun which is in a circular orbit
         around the Galactic centre.  The Sun\texttt{'} s motion with respect to
         the dynamical LSR is called the \texttt{"} peculiar\texttt{"}  solar motion.  To
         obtain a radial velocity correction with respect to the
         dynamical LSR use the routine palRvlsrd.
      }
   }
   \sstdiytopic{
      Reference
   }{
      \sstitemlist{

         \sstitem
          Delhaye (1965), in \texttt{"} Stars and Stellar Systems\texttt{"} , vol 5, p73.
      }
   }
}
\sstroutine{
   palSubet
}{
   Remove the E-terms from a pre IAU 1976 catalogue RA,Dec
}{
   \sstdescription{
      Remove the E-terms (elliptic component of annual aberration)
      from a pre IAU 1976 catalogue RA,Dec to give a mean place.
   }
   \sstinvocation{
      void palSubet ( double rc, double dc, double eq,
                      double $*$rm, double $*$dm );
   }
   \sstarguments{
      \sstsubsection{
         rc = double (Given)
      }{
         RA with E-terms included (radians)
      }
      \sstsubsection{
         dc = double (Given)
      }{
         Dec with E-terms included (radians)
      }
      \sstsubsection{
         eq = double (Given)
      }{
         Besselian epoch of mean equator and equinox
      }
      \sstsubsection{
         rm = double $*$ (Returned)
      }{
         RA without E-terms (radians)
      }
      \sstsubsection{
         dm = double $*$ (Returned)
      }{
         Dec without E-terms (radians)
      }
   }
   \sstnotes{
      Most star positions from pre-1984 optical catalogues (or
      derived from astrometry using such stars) embody the
      E-terms.  This routine converts such a position to a
      formal mean place (allowing, for example, comparison with a
      pulsar timing position).
   }
   \sstdiytopic{
      See Also
   }{
      Explanatory Supplement to the Astronomical Ephemeris,
      section 2D, page 48.
   }
}
\sstroutine{
   palSupgal
}{
   Convert from supergalactic to galactic coordinates
}{
   \sstdescription{
      Transformation from de Vaucouleurs supergalactic coordinates
      to IAU 1958 galactic coordinates
   }
   \sstinvocation{
      void palSupgal ( double dsl, double dsb, double $*$dl, double $*$db );
   }
   \sstarguments{
      \sstsubsection{
         dsl = double (Given)
      }{
         Supergalactic longitude.
      }
      \sstsubsection{
         dsb = double (Given)
      }{
         Supergalactic latitude.
      }
      \sstsubsection{
         dl = double $*$ (Returned)
      }{
         Galactic longitude.
      }
      \sstsubsection{
         db = double $*$ (Returned)
      }{
         Galactic latitude.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          de Vaucouleurs, de Vaucouleurs, \& Corwin, Second Reference
           Catalogue of Bright Galaxies, U. Texas, page 8.

         \sstitem
          Systems \& Applied Sciences Corp., Documentation for the
           machine-readable version of the above catalogue,
           Contract NAS 5-26490.

      }
      (These two references give different values for the galactic
      longitude of the supergalactic origin.  Both are wrong;  the
      correct value is L2=137.37.)
   }
}
\sstroutine{
   palUe2el
}{
   Universal elements to heliocentric osculating elements
}{
   \sstdescription{
      Transform universal elements into conventional heliocentric
      osculating elements.
   }
   \sstinvocation{
      void palUe2el ( const double u[13], int jformr,
                      int $*$jform, double $*$epoch, double $*$orbinc,
                      double $*$anode, double $*$perih, double $*$aorq, double $*$e,
                      double $*$aorl, double $*$dm, int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         u = const double [13] (Given)
      }{
         Universal orbital elements (Note 1)
             (0)  combined mass (M$+$m)
             (1)  total energy of the orbit (alpha)
             (2)  reference (osculating) epoch (t0)
           (3-5)  position at reference epoch (r0)
           (6-8)  velocity at reference epoch (v0)
             (9)  heliocentric distance at reference epoch
            (10)  r0.v0
            (11)  date (t)
            (12)  universal eccentric anomaly (psi) of date, approx
      }
      \sstsubsection{
         jformr = int (Given)
      }{
         Requested element set (1-3; Note 3)
      }
      \sstsubsection{
         jform = int $*$ (Returned)
      }{
         Element set actually returned (1-3; Note 4)
      }
      \sstsubsection{
         epoch = double $*$ (Returned)
      }{
         Epoch of elements (TT MJD)
      }
      \sstsubsection{
         orbinc = double $*$ (Returned)
      }{
         inclination (radians)
      }
      \sstsubsection{
         anode = double $*$ (Returned)
      }{
         longitude of the ascending node (radians)
      }
      \sstsubsection{
         perih = double $*$ (Returned)
      }{
         longitude or argument of perihelion (radians)
      }
      \sstsubsection{
         aorq = double $*$ (Returned)
      }{
         mean distance or perihelion distance (AU)
      }
      \sstsubsection{
         e = double $*$ (Returned)
      }{
         eccentricity
      }
      \sstsubsection{
         aorl = double $*$ (Returned)
      }{
         mean anomaly or longitude (radians, JFORM=1,2 only)
      }
      \sstsubsection{
         dm = double $*$ (Returned)
      }{
         daily motion (radians, JFORM=1 only)
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:  0 = OK
         \sstitemlist{

            \sstitem
                     1 = illegal combined mass

            \sstitem
                     2 = illegal JFORMR

            \sstitem
                     3 = position/velocity out of range
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The \texttt{"} universal\texttt{"}  elements are those which define the orbit for the
           purposes of the method of universal variables (see reference 2).
           They consist of the combined mass of the two bodies, an epoch,
           and the position and velocity vectors (arbitrary reference frame)
           at that epoch.  The parameter set used here includes also various
           quantities that can, in fact, be derived from the other
           information.  This approach is taken to avoiding unnecessary
           computation and loss of accuracy.  The supplementary quantities
           are (i) alpha, which is proportional to the total energy of the
           orbit, (ii) the heliocentric distance at epoch, (iii) the
           outwards component of the velocity at the given epoch, (iv) an
           estimate of psi, the \texttt{"} universal eccentric anomaly\texttt{"}  at a given
           date and (v) that date.

         \sstitem
          The universal elements are with respect to the mean equator and
           equinox of epoch J2000.  The orbital elements produced are with
           respect to the J2000 ecliptic and mean equinox.

         \sstitem
          Three different element-format options are supported:

      }
         Option JFORM=1, suitable for the major planets:

         EPOCH  = epoch of elements (TT MJD)
         ORBINC = inclination i (radians)
         ANODE  = longitude of the ascending node, big omega (radians)
         PERIH  = longitude of perihelion, curly pi (radians)
         AORQ   = mean distance, a (AU)
         E      = eccentricity, e
         AORL   = mean longitude L (radians)
         DM     = daily motion (radians)

         Option JFORM=2, suitable for minor planets:

         EPOCH  = epoch of elements (TT MJD)
         ORBINC = inclination i (radians)
         ANODE  = longitude of the ascending node, big omega (radians)
         PERIH  = argument of perihelion, little omega (radians)
         AORQ   = mean distance, a (AU)
         E      = eccentricity, e
         AORL   = mean anomaly M (radians)

         Option JFORM=3, suitable for comets:

         EPOCH  = epoch of perihelion (TT MJD)
         ORBINC = inclination i (radians)
         ANODE  = longitude of the ascending node, big omega (radians)
         PERIH  = argument of perihelion, little omega (radians)
         AORQ   = perihelion distance, q (AU)
         E      = eccentricity, e

      \sstitemlist{

         \sstitem
          It may not be possible to generate elements in the form
           requested through JFORMR.  The caller is notified of the form
           of elements actually returned by means of the JFORM argument:

      }
         JFORMR   JFORM     meaning

           1        1       OK - elements are in the requested format
           1        2       never happens
           1        3       orbit not elliptical

           2        1       never happens
           2        2       OK - elements are in the requested format
           2        3       orbit not elliptical

           3        1       never happens
           3        2       never happens
           3        3       OK - elements are in the requested format

      \sstitemlist{

         \sstitem
          The arguments returned for each value of JFORM (cf Note 6: JFORM
           may not be the same as JFORMR) are as follows:

      }
          JFORM         1              2              3
          EPOCH         t0             t0             T
          ORBINC        i              i              i
          ANODE         Omega          Omega          Omega
          PERIH         curly pi       omega          omega
          AORQ          a              a              q
          E             e              e              e
          AORL          L              M              -
          DM            n              -              -

      where:

          t0           is the epoch of the elements (MJD, TT)
          T              \texttt{"}     epoch of perihelion (MJD, TT)
          i              \texttt{"}     inclination (radians)
          Omega          \texttt{"}     longitude of the ascending node (radians)
          curly pi       \texttt{"}     longitude of perihelion (radians)
          omega          \texttt{"}     argument of perihelion (radians)
          a              \texttt{"}     mean distance (AU)
          q              \texttt{"}     perihelion distance (AU)
          e              \texttt{"}     eccentricity
          L              \texttt{"}     longitude (radians, 0-2pi)
          M              \texttt{"}     mean anomaly (radians, 0-2pi)
          n              \texttt{"}     daily motion (radians)
      \sstitemlist{

         \sstitem
              means no value is set

         \sstitem
          At very small inclinations, the longitude of the ascending node
           ANODE becomes indeterminate and under some circumstances may be
           set arbitrarily to zero.  Similarly, if the orbit is close to
           circular, the true anomaly becomes indeterminate and under some
           circumstances may be set arbitrarily to zero.  In such cases,
           the other elements are automatically adjusted to compensate,
           and so the elements remain a valid description of the orbit.
      }
   }
   \sstdiytopic{
      See Also
   }{
      \sstitemlist{

         \sstitem
          Sterne, Theodore E., \texttt{"} An Introduction to Celestial Mechanics\texttt{"} ,
           Interscience Publishers Inc., 1960.  Section 6.7, p199.

         \sstitem
          Everhart, E. \& Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
      }
   }
}
\sstroutine{
   palUe2pv
}{
   Heliocentric position and velocity of a planet, asteroid or comet, from universal elements
}{
   \sstdescription{
      Heliocentric position and velocity of a planet, asteroid or comet,
      starting from orbital elements in the \texttt{"} universal variables\texttt{"}  form.
   }
   \sstinvocation{
      void palUe2pv( double date, double u[13], double pv[6], int $*$jstat );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT Modified Julian date (JD-2400000.5).
      }
      \sstsubsection{
         u = double [13] (Given \& Returned)
      }{
          Universal orbital elements (updated, see note 1)
          given    (0)   combined mass (M$+$m)
            \texttt{"}       (1)   total energy of the orbit (alpha)
            \texttt{"}       (2)   reference (osculating) epoch (t0)
            \texttt{"}     (3-5)   position at reference epoch (r0)
            \texttt{"}     (6-8)   velocity at reference epoch (v0)
            \texttt{"}       (9)   heliocentric distance at reference epoch
            \texttt{"}      (10)   r0.v0
         returned (11)   date (t)
            \texttt{"}      (12)   universal eccentric anomaly (psi) of date
      }
      \sstsubsection{
         pv = double [6] (Returned)
      }{
         Position (AU) and velocity (AU/s)
      }
      \sstsubsection{
         jstat = int $*$ (Returned)
      }{
         status:  0 = OK
         \sstitemlist{

            \sstitem
                     1 = radius vector zero

            \sstitem
                     2 = failed to converge
         }
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The \texttt{"} universal\texttt{"}  elements are those which define the orbit for the
           purposes of the method of universal variables (see reference).
           They consist of the combined mass of the two bodies, an epoch,
           and the position and velocity vectors (arbitrary reference frame)
           at that epoch.  The parameter set used here includes also various
           quantities that can, in fact, be derived from the other
           information.  This approach is taken to avoiding unnecessary
           computation and loss of accuracy.  The supplementary quantities
           are (i) alpha, which is proportional to the total energy of the
           orbit, (ii) the heliocentric distance at epoch, (iii) the
           outwards component of the velocity at the given epoch, (iv) an
           estimate of psi, the \texttt{"} universal eccentric anomaly\texttt{"}  at a given
           date and (v) that date.

         \sstitem
          The companion routine is palEl2ue.  This takes the conventional
           orbital elements and transforms them into the set of numbers
           needed by the present routine.  A single prediction requires one
           one call to palEl2ue followed by one call to the present routine;
           for convenience, the two calls are packaged as the routine
           palPlanel.  Multiple predictions may be made by again
           calling palEl2ue once, but then calling the present routine
           multiple times, which is faster than multiple calls to palPlanel.

         \sstitem
          It is not obligatory to use palEl2ue to obtain the parameters.
           However, it should be noted that because palEl2ue performs its
           own validation, no checks on the contents of the array U are made
           by the present routine.

         \sstitem
          DATE is the instant for which the prediction is required.  It is
           in the TT timescale (formerly Ephemeris Time, ET) and is a
           Modified Julian Date (JD-2400000.5).

         \sstitem
          The universal elements supplied in the array U are in canonical
           units (solar masses, AU and canonical days).  The position and
           velocity are not sensitive to the choice of reference frame.  The
           palEl2ue routine in fact produces coordinates with respect to the
           J2000 equator and equinox.

         \sstitem
          The algorithm was originally adapted from the EPHSLA program of
           D.H.P.Jones (private communication, 1996).  The method is based
           on Stumpff\texttt{'} s Universal Variables.

         \sstitem
          Reference:  Everhart, E. \& Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
      }
   }
}
\sstroutine{
   palUnpcd
}{
   Remove pincushion/barrel distortion
}{
   \sstdescription{
      Remove pincushion/barrel distortion from a distorted [x,y] to give
      tangent-plane [x,y].
   }
   \sstinvocation{
      palUnpcd( double disco, double $*$ x, double $*$ y );
   }
   \sstarguments{
      \sstsubsection{
         disco = double (Given)
      }{
         Pincushion/barrel distortion coefficient.
      }
      \sstsubsection{
         x = double $*$ (Given \& Returned)
      }{
         On input the distorted X coordinate, on output
         the tangent-plane X coordinate.
      }
      \sstsubsection{
         y = double $*$ (Given \& Returned)
      }{
         On input the distorted Y coordinate, on output
         the tangent-plane Y coordinate.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The distortion is of the form RP = R$*$(1$+$C$*$R$\wedge$2), where R is
           the radial distance from the tangent point, C is the DISCO
           argument, and RP is the radial distance in the presence of
           the distortion.

         \sstitem
          For pincushion distortion, C is $+$ve;  for barrel distortion,
           C is -ve.

         \sstitem
          For X,Y in \texttt{"} radians\texttt{"}  - units of one projection radius,
           which in the case of a photograph is the focal length of
           the camera - the following DISCO values apply:

      }
            Geometry          DISCO

            astrograph         0.0
            Schmidt           -0.3333
            AAT PF doublet  $+$147.069
            AAT PF triplet  $+$178.585
            AAT f/8          $+$21.20
            JKT f/8          $+$13.32

      \sstitemlist{

         \sstitem
          The present routine is a rigorous inverse of the companion
           routine palPcd.  The expression for RP in Note 1 is rewritten
           in the form x$\wedge$3$+$a$*$x$+$b=0 and solved by standard techniques.

         \sstitem
          Cases where the cubic has multiple real roots can sometimes
           occur, corresponding to extreme instances of barrel distortion
           where up to three different undistorted [X,Y]s all produce the
           same distorted [X,Y].  However, only one solution is returned,
           the one that produces the smallest change in [X,Y].
      }
   }
   \sstdiytopic{
      See Also
   }{
      palPcd
   }
}
\sstroutine{
   palVers
}{
   Obtain PAL version number
}{
   \sstdescription{
      Retrieve the PAL version number as a string in the form \texttt{"} A.B.C\texttt{"}
      and as an integer (major$*$1e6$+$minor$*$1e3$+$release).
   }
   \sstinvocation{
      int palVers( char $*$verstring, size\_t verlen );
   }
   \sstarguments{
      \sstsubsection{
         verstring = char $*$ (Returned)
      }{
         Buffer to receive version string of the form \texttt{"} A.B.C\texttt{"} . Can be NULL.
      }
      \sstsubsection{
         verlen = size\_t (Given)
      }{
         Allocated size of \texttt{"} verstring\texttt{"}  including nul. Version string
         will be truncated if it does not fit in buffer.
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         vernum = int (Returned)
      }{
         Version number as an integer.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Note that this API does not match the slaVers API.
      }
   }
}
\sstroutine{
   palXy2xy
}{
   Transform one [x,y] into another using a linear model of the type
   produced by the palFitxy routine
}{
   \sstdescription{
      The model relates two sets of [x,y] coordinates as follows.
      Naming the elements of coeffs:
      \newline
      \newline
      \hspace*{1.5 em}
         coeffs[0] = A
      \newline
      \hspace*{1.5 em}
         coeffs[1] = B
      \newline
      \hspace*{1.5 em}
         coeffs[2] = C
      \newline
      \hspace*{1.5 em}
         coeffs[3] = D
      \newline
      \hspace*{1.5 em}
         coeffs[4] = E
      \newline
      \hspace*{1.5 em}
         coeffs[5] = F
      \newline
      \newline
      the present routine performs the transformation:
      \newline
      \newline
      \hspace*{1.5 em}
         x2 = A $+$ B $*$ x1 $+$ C $*$ y1
      \newline
      \hspace*{1.5 em}
         y2 = D $+$ E $*$ x1 $+$ F $*$ y1
      \newline
      \newline
   }
   \sstinvocation{
      palXy2xy ( double x1, double y1, double coeffs[6],
                 double $*$x2, double $*$y2)
   }
   \sstarguments{
      \sstsubsection{
         x1 = double (Given)
      }{
         x-coordinate
      }
      \sstsubsection{
         y1 = double (Given)
      }{
         y-coordinate
      }
      \sstsubsection{
         coeffs = double[6] (Given)
      }{
         transformation coefficients (see note)
      }
      \sstsubsection{
         x2 = double (Returned)
      }{
         x-coordinate
      }
      \sstsubsection{
         y2 = double (Returned)
      }{
         y-coordinate
      }
   }
   \sstdiytopic{
      See also
   }{
      palFitxy, palPxy, palInvf and palDcmpf
   }
}
\sstroutine{
   palZd
}{
   HA, Dec to Zenith Distance
}{
   \sstdescription{
      Calculate the zenith distance for a given hour angle and declination
      using the SOFA/ERFA routine eraHd2ae.
   }
   \sstinvocation{
      double palZd( double ha, double dec, double phi );
   }
   \sstarguments{
      \sstsubsection{
         ha = double (Given)
      }{
         Hour angle (radians)
      }
      \sstsubsection{
         dec = double (Given)
      }{
         Declination (radians)
      }
      \sstsubsection{
         phi = double (Given)
      }{
         Observatory latitude (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraHd2ae and converts the returned elevation to zenith
           distance.  If elevation itself, or azimuth, is desired then
           that routine could instead be called directly.
      }
   }
}
\sstroutine{
   palCldj
}{
   Gregorian Calendar to Modified Julian Date
}{
   \sstdescription{
      Gregorian calendar to Modified Julian Date.
   }
   \sstinvocation{
      palCldj( int iy, int im, int id, double $*$djm, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         iy = int (Given)
      }{
         Year in Gregorian calendar
      }
      \sstsubsection{
         im = int (Given)
      }{
         Month in Gregorian calendar
      }
      \sstsubsection{
         id = int (Given)
      }{
         Day in Gregorian calendar
      }
      \sstsubsection{
         djm = double $*$ (Returned)
      }{
         Modified Julian Date (JD-2400000.5) for 0 hrs
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         status: 0 = OK, 1 = bad year (MJD not computed),
         2 = bad month (MJD not computed), 3 = bad day (MJD computed).
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraCal2jd(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDbear
}{
   Bearing (position angle) of one point on a sphere relative to another
}{
   \sstdescription{
      Bearing (position angle) of one point in a sphere relative to another.
   }
   \sstinvocation{
      pa = palDbear( double a1, double b1, double a2, double b2 );
   }
   \sstarguments{
      \sstsubsection{
         a1 = double (Given)
      }{
         Longitude of point 1 (e.g. RA) in radians.
      }
      \sstsubsection{
         b1 = double (Given)
      }{
         Latitude of point 1 (e.g. Dec) in radians.
      }
      \sstsubsection{
         a2 = double (Given)
      }{
         Longitude of point 2 in radians.
      }
      \sstsubsection{
         b2 = double (Given)
      }{
         Latitude of point 2 in radians.
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         The result is the bearing (position angle), in radians, of point
      }{
      }
      \sstsubsection{
         A2,B2 as seen from point A1,B1.  It is in the range $+$/- pi.  If
      }{
      }
      \sstsubsection{
         A2,B2 is due east of A1,B1 the bearing is $+$pi/2.  Zero is returned
      }{
      }
      \sstsubsection{
         if the two points are coincident.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraPas(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDaf2r
}{
   Convert degrees, arcminutes, arcseconds to radians
}{
   \sstdescription{
      Convert degrees, arcminutes, arcseconds to radians.
   }
   \sstinvocation{
      palDaf2r( int ideg, int iamin, double asec, double $*$rad, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         ideg = int (Given)
      }{
         Degrees.
      }
      \sstsubsection{
         iamin = int (Given)
      }{
         Arcminutes.
      }
      \sstsubsection{
         iasec = double (Given)
      }{
         Arcseconds.
      }
      \sstsubsection{
         rad = double $*$ (Returned)
      }{
         Angle in radians.
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         Status: 0 = OK, 1 = \texttt{"} ideg\texttt{"}  out of range 0-359,
                 2 = \texttt{"} iamin\texttt{"}  outside of range 0-59,
                 2 = \texttt{"} asec\texttt{"}  outside range 0-59.99999
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraAf2a(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDav2m
}{
   Form the rotation matrix corresponding to a given axial vector
}{
   \sstdescription{
      A rotation matrix describes a rotation about some arbitrary axis,
      called the Euler axis.  The \texttt{"} axial vector\texttt{"}  supplied to this routine
      has the same direction as the Euler axis, and its magnitude is the
      amount of rotation in radians.
   }
   \sstinvocation{
      palDav2m( double axvec[3], double rmat[3][3] );
   }
   \sstarguments{
      \sstsubsection{
         axvec = double [3] (Given)
      }{
         Axial vector (radians)
      }
      \sstsubsection{
         rmat = double [3][3] (Returned)
      }{
         Rotation matrix.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraRv2m(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDcc2s
}{
   Cartesian to spherical coordinates
}{
   \sstdescription{
      The spherical coordinates are longitude ($+$ve anticlockwise looking
      from the $+$ve latitude pole) and latitude.  The Cartesian coordinates
      are right handed, with the x axis at zero longitude and latitude, and
      the z axis at the $+$ve latitude pole.
   }
   \sstinvocation{
      palDcc2s( double v[3], double $*$a, double $*$b );
   }
   \sstarguments{
      \sstsubsection{
         v = double [3] (Given)
      }{
         x, y, z vector.
      }
      \sstsubsection{
         a = double $*$ (Returned)
      }{
         Spherical coordinate (radians)
      }
      \sstsubsection{
         b = double $*$ (Returned)
      }{
         Spherical coordinate (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraC2s(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDcs2c
}{
   Spherical coordinates to direction cosines
}{
   \sstdescription{
      The spherical coordinates are longitude ($+$ve anticlockwise looking
      from the $+$ve latitude pole) and latitude.  The Cartesian coordinates
      are right handed, with the x axis at zero longitude and latitude, and
      the z axis at the $+$ve latitude pole.
   }
   \sstinvocation{
      palDcs2c( double a, double b, double v[3] );
   }
   \sstarguments{
      \sstsubsection{
         a = double (Given)
      }{
         Spherical coordinate in radians (ra, long etc).
      }
      \sstsubsection{
         b = double (Given)
      }{
         Spherical coordinate in radians (dec, lat etc).
      }
      \sstsubsection{
         v = double [3] (Returned)
      }{
         x, y, z vector
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraS2c(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDd2tf
}{
   Convert an interval in days into hours, minutes, seconds
}{
   \sstdescription{
      Convert and interval in days into hours, minutes, seconds.
   }
   \sstinvocation{
      palDd2tf( int ndp, double days, char $*$sign, int ihmsf[4] );
   }
   \sstarguments{
      \sstsubsection{
         ndp = int (Given)
      }{
         Number of decimal places of seconds
      }
      \sstsubsection{
         days = double (Given)
      }{
         Interval in days
      }
      \sstsubsection{
         sign = char $*$ (Returned)
      }{
         \texttt{'} $+$\texttt{'}  or \texttt{'} -\texttt{'}  (single character, not string)
      }
      \sstsubsection{
         ihmsf = int [4] (Returned)
      }{
         Hours, minutes, seconds, fraction
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraD2tf(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDimxv
}{
   Perform the 3-D backward unitary transformation
}{
   \sstdescription{
      Perform the 3-D backward unitary transformation.
   }
   \sstinvocation{
      palDimxv( double dm[3][3], double va[3], double vb[3] );
   }
   \sstarguments{
      \sstsubsection{
         dm = double [3][3] (Given)
      }{
         Matrix
      }
      \sstsubsection{
         va = double [3] (Given)
      }{
         vector
      }
      \sstsubsection{
         vb = double [3] (Returned)
      }{
         Result vector
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraTrxp(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDm2av
}{
   From a rotation matrix, determine the corresponding axial vector
}{
   \sstdescription{
      A rotation matrix describes a rotation about some arbitrary axis,
      called the Euler axis.  The \texttt{"} axial vector\texttt{"}  returned by this routine
      has the same direction as the Euler axis, and its magnitude is the
      amount of rotation in radians.  (The magnitude and direction can be
      separated by means of the routine palDvn.)
   }
   \sstinvocation{
      palDm2av( double rmat[3][3], double axvec[3] );
   }
   \sstarguments{
      \sstsubsection{
         rmat = double [3][3] (Given)
      }{
         Rotation matrix
      }
      \sstsubsection{
         axvec = double [3] (Returned)
      }{
         Axial vector (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraRm2v(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDjcl
}{
   Modified Julian Date to Gregorian year, month, day and fraction of day
}{
   \sstdescription{
      Modified Julian Date to Gregorian year, month, day and fraction of day.
   }
   \sstinvocation{
      palDjcl( double djm, int $*$iy, int $*$im, int $*$id, double $*$fd, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         djm = double (Given)
      }{
         modified Julian Date (JD-2400000.5)
      }
      \sstsubsection{
         iy = int $*$ (Returned)
      }{
         year
      }
      \sstsubsection{
         im = int $*$ (Returned)
      }{
         month
      }
      \sstsubsection{
         id = int $*$ (Returned)
      }{
         day
      }
      \sstsubsection{
         fd = double $*$ (Returned)
      }{
         Fraction of day.
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraJd2cal(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDmxm
}{
   Product of two 3x3 matrices
}{
   \sstdescription{
      Product of two 3x3 matrices.
   }
   \sstinvocation{
      palDmxm( double a[3][3], double b[3][3], double c[3][3] );
   }
   \sstarguments{
      \sstsubsection{
         a = double [3][3] (Given)
      }{
         Matrix
      }
      \sstsubsection{
         b = double [3][3] (Given)
      }{
         Matrix
      }
      \sstsubsection{
         c = double [3][3] (Returned)
      }{
         Matrix result
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraRxr(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDmxv
}{
   Performs the 3-D forward unitary transformation
}{
   \sstdescription{
      Performs the 3-D forward unitary transformation.
   }
   \sstinvocation{
      palDmxv( double dm[3][3], double va[3], double vb[3] );
   }
   \sstarguments{
      \sstsubsection{
         dm = double [3][3] (Given)
      }{
         matrix
      }
      \sstsubsection{
         va = double [3] (Given)
      }{
         vector
      }
      \sstsubsection{
         dp = double [3] (Returned)
      }{
         result vector
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraRxp(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDpav
}{
   Position angle of one celestial direction with respect to another
}{
   \sstdescription{
      Position angle of one celestial direction with respect to another.
   }
   \sstinvocation{
      pa = palDpav( double v1[3], double v2[3] );
   }
   \sstarguments{
      \sstsubsection{
         v1 = double [3] (Given)
      }{
         direction cosines of one point.
      }
      \sstsubsection{
         v2 = double [3] (Given)
      }{
         direction cosines of the other point.
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         The result is the bearing (position angle), in radians, of point
      }{
      }
      \sstsubsection{
         V2 with respect to point V1.  It is in the range $+$/- pi.  The
      }{
      }
      \sstsubsection{
         sense is such that if V2 is a small distance east of V1, the
      }{
      }
      \sstsubsection{
         bearing is about $+$pi/2.  Zero is returned if the two points
      }{
      }
      \sstsubsection{
         are coincident.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The coordinate frames correspond to RA,Dec, Long,Lat etc.

         \sstitem
          Uses eraPap(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDr2af
}{
   Convert an angle in radians to degrees, arcminutes, arcseconds
}{
   \sstdescription{
      Convert an angle in radians to degrees, arcminutes, arcseconds.
   }
   \sstinvocation{
      palDr2af( int ndp, double angle, char $*$sign, int idmsf[4] );
   }
   \sstarguments{
      \sstsubsection{
         ndp = int (Given)
      }{
         number of decimal places of arcseconds
      }
      \sstsubsection{
         angle = double (Given)
      }{
         angle in radians
      }
      \sstsubsection{
         sign = char $*$ (Returned)
      }{
         \texttt{'} $+$\texttt{'}  or \texttt{'} -\texttt{'}  (single character)
      }
      \sstsubsection{
         idmsf = int [4] (Returned)
      }{
         Degrees, arcminutes, arcseconds, fraction
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraA2af(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDr2tf
}{
   Convert an angle in radians to hours, minutes, seconds
}{
   \sstdescription{
      Convert an angle in radians to hours, minutes, seconds.
   }
   \sstinvocation{
      palDr2tf ( int ndp, double angle, char $*$sign, int ihmsf[4] );
   }
   \sstarguments{
      \sstsubsection{
         ndp = int (Given)
      }{
         number of decimal places of arcseconds
      }
      \sstsubsection{
         angle = double (Given)
      }{
         angle in radians
      }
      \sstsubsection{
         sign = char $*$ (Returned)
      }{
         \texttt{'} $+$\texttt{'}  or \texttt{'} -\texttt{'}  (single character)
      }
      \sstsubsection{
         idmsf = int [4] (Returned)
      }{
         Hours, minutes, seconds, fraction
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraA2tf(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDranrm
}{
   Normalize angle into range 0-2 pi
}{
   \sstdescription{
      Normalize angle into range 0-2 pi.
   }
   \sstinvocation{
      norm = palDranrm( double angle );
   }
   \sstarguments{
      \sstsubsection{
         angle = double (Given)
      }{
         angle in radians
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Angle expressed in the range 0-2 pi
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraAnp(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDsep
}{
   Angle between two points on a sphere
}{
   \sstdescription{
      Angle between two points on a sphere.
   }
   \sstinvocation{
      ang = palDsep( double a1, double b1, double a2, double b2 );
   }
   \sstarguments{
      \sstsubsection{
         a1 = double (Given)
      }{
         Spherical coordinate of one point (radians)
      }
      \sstsubsection{
         b1 = double (Given)
      }{
         Spherical coordinate of one point (radians)
      }
      \sstsubsection{
         a2 = double (Given)
      }{
         Spherical coordinate of other point (radians)
      }
      \sstsubsection{
         b2 = double (Given)
      }{
         Spherical coordinate of other point (radians)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Angle, in radians, between the two points. Always positive.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          The spherical coordinates are [RA,Dec], [Long,Lat] etc, in radians.

         \sstitem
          Uses eraSeps(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDsepv
}{
   Angle between two vectors
}{
   \sstdescription{
      Angle between two vectors.
   }
   \sstinvocation{
      ang = palDsepv( double v1[3], double v2[3] );
   }
   \sstarguments{
      \sstsubsection{
         v1 = double [3] (Given)
      }{
         First vector
      }
      \sstsubsection{
         v2 = double [3] (Given)
      }{
         Second vector
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Angle, in radians, between the two points. Always positive.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraSepp(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDtf2d
}{
   Convert hours, minutes, seconds to days
}{
   \sstdescription{
      Convert hours, minutes, seconds to days.
   }
   \sstinvocation{
      palDtf2d( int ihour, int imin, double sec, double $*$days, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         ihour = int (Given)
      }{
         Hours
      }
      \sstsubsection{
         imin = int (Given)
      }{
         Minutes
      }
      \sstsubsection{
         sec = double (Given)
      }{
         Seconds
      }
      \sstsubsection{
         days = double $*$ (Returned)
      }{
         Interval in days
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         status: 0 = ok, 1 = ihour outside range 0-23,
         2 = imin outside range 0-59, 3 = sec outside range 0-59.999...
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraTf2d(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDtf2r
}{
   Convert hours, minutes, seconds to radians
}{
   \sstdescription{
      Convert hours, minutes, seconds to radians.
   }
   \sstinvocation{
      palDtf2r( int ihour, int imin, double sec, double $*$rad, int $*$j );
   }
   \sstarguments{
      \sstsubsection{
         ihour = int (Given)
      }{
         Hours
      }
      \sstsubsection{
         imin = int (Given)
      }{
         Minutes
      }
      \sstsubsection{
         sec = double (Given)
      }{
         Seconds
      }
      \sstsubsection{
         days = double $*$ (Returned)
      }{
         Angle in radians
      }
      \sstsubsection{
         j = int $*$ (Returned)
      }{
         status: 0 = ok, 1 = ihour outside range 0-23,
         2 = imin outside range 0-59, 3 = sec outside range 0-59.999...
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraTf2a(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDvdv
}{
   Scalar product of two 3-vectors
}{
   \sstdescription{
      Scalar product of two 3-vectors.
   }
   \sstinvocation{
      prod = palDvdv ( double va[3], double vb[3] );
   }
   \sstarguments{
      \sstsubsection{
         va = double [3] (Given)
      }{
         First vector
      }
      \sstsubsection{
         vb = double [3] (Given)
      }{
         Second vector
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Scalar product va.vb
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraPdp(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDvn
}{
   Normalizes a 3-vector also giving the modulus
}{
   \sstdescription{
      Normalizes a 3-vector also giving the modulus.
   }
   \sstinvocation{
      palDvn( double v[3], double uv[3], double $*$vm );
   }
   \sstarguments{
      \sstsubsection{
         v = double [3] (Given)
      }{
         vector
      }
      \sstsubsection{
         uv = double [3] (Returned)
      }{
         unit vector in direction of \texttt{"} v\texttt{"}
      }
      \sstsubsection{
         vm = double $*$ (Returned)
      }{
         modulus of \texttt{"} v\texttt{"}
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraPn(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palDvxv
}{
   Vector product of two 3-vectors
}{
   \sstdescription{
      Vector product of two 3-vectors.
   }
   \sstinvocation{
      palDvxv( double va[3], double vb[3], double vc[3] );
   }
   \sstarguments{
      \sstsubsection{
         va = double [3] (Given)
      }{
         First vector
      }
      \sstsubsection{
         vb = double [3] (Given)
      }{
         Second vector
      }
      \sstsubsection{
         vc = double [3] (Returned)
      }{
         Result vector
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraPxp(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palEpb
}{
   Conversion of modified Julian Data to Besselian Epoch
}{
   \sstdescription{
      Conversion of modified Julian Data to Besselian Epoch.
   }
   \sstinvocation{
      epb = palEpb ( double date );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         Modified Julian Date (JD - 2400000.5)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Besselian epoch.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraEpb(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palEpb2d
}{
   Conversion of Besselian Epoch to Modified Julian Date
}{
   \sstdescription{
      Conversion of Besselian Epoch to Modified Julian Date.
   }
   \sstinvocation{
      mjd = palEpb2d ( double epb );
   }
   \sstarguments{
      \sstsubsection{
         epb = double (Given)
      }{
         Besselian Epoch
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Modified Julian Date (JD - 2400000.5)
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraEpb2jd(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palEpj
}{
   Conversion of Modified Julian Date to Julian Epoch
}{
   \sstdescription{
      Conversion of Modified Julian Date to Julian Epoch.
   }
   \sstinvocation{
      epj = palEpj ( double date );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         Modified Julian Date (JD - 2400000.5)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         The Julian Epoch.
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraEpj(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palEpj2d
}{
   Conversion of Julian Epoch to Modified Julian Date
}{
   \sstdescription{
      Conversion of Julian Epoch to Modified Julian Date.
   }
   \sstinvocation{
      mjd = palEpj2d ( double epj );
   }
   \sstarguments{
      \sstsubsection{
         epj = double (Given)
      }{
         Julian Epoch.
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Modified Julian Date (JD - 2400000.5)
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraEpj2d(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palEqeqx
}{
   Equation of the equinoxes (IAU 2000/2006)
}{
   \sstdescription{
      Equation of the equinoxes (IAU 2000/2006).
   }
   \sstinvocation{
      palEqeqx( double date );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         TT as Modified Julian Date (JD-400000.5)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraEe06a(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palFk5hz
}{
   Transform an FK5 (J2000) star position into the frame of the
   Hipparcos catalogue
}{
   \sstdescription{
      Transform an FK5 (J2000) star position into the frame of the
      Hipparcos catalogue.
   }
   \sstinvocation{
      palFk5hz ( double r5, double d5, double epoch,
                 double $*$rh, double $*$dh );
   }
   \sstarguments{
      \sstsubsection{
         r5 = double (Given)
      }{
         FK5 RA (radians), equinox J2000, epoch \texttt{"} epoch\texttt{"}
      }
      \sstsubsection{
         d5 = double (Given)
      }{
         FK5 dec (radians), equinox J2000, epoch \texttt{"} epoch\texttt{"}
      }
      \sstsubsection{
         epoch = double (Given)
      }{
         Julian epoch
      }
      \sstsubsection{
         rh = double $*$ (Returned)
      }{
         RA (radians)
      }
      \sstsubsection{
         dh = double $*$ (Returned)
      }{
         Dec (radians)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Assumes zero Hipparcos proper motion.

         \sstitem
          Uses eraEpj2jd() and eraFk5hz.
           See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palGmst
}{
   Greenwich mean sidereal time (consistent with IAU 2006 precession)
}{
   \sstdescription{
      Greenwich mean sidereal time (consistent with IAU 2006 precession).
   }
   \sstinvocation{
      mst = palGmst ( double ut1 );
   }
   \sstarguments{
      \sstsubsection{
         ut1 = double (Given)
      }{
         Universal time (UT1) expressed as modified Julian Date (JD-2400000.5)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Greenwich mean sidereal time
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraGmst06(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palGmsta
}{
   Greenwich mean sidereal time (consistent with IAU 2006 precession)
}{
   \sstdescription{
      Greenwich mean sidereal time (consistent with IAU 2006 precession).
   }
   \sstinvocation{
      mst = palGmsta ( double date, double ut1 );
   }
   \sstarguments{
      \sstsubsection{
         date = double (Given)
      }{
         UT1 date (MJD: integer part of JD-2400000.5)
      }
      \sstsubsection{
         ut1 = double (Given)
      }{
         UT1 time (fraction of a day)
      }
   }
   \sstreturnedvalue{
      \sstsubsection{
         Greenwich mean sidereal time (in range 0 to 2 pi)
      }{
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          For best accuracy use eraGmst06() directly.

         \sstitem
          Uses eraGmst06(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palHfk5z
}{
   Hipparcos star position to FK5 J2000
}{
   \sstdescription{
      Transform a Hipparcos star position into FK5 J2000, assuming
      zero Hipparcos proper motion.
   }
   \sstinvocation{
      palHfk5z( double rh, double dh, double epoch,
                double $*$r5, double $*$d5, double $*$dr5, double $*$dd5 );
   }
   \sstarguments{
      \sstsubsection{
         rh = double (Given)
      }{
         Hipparcos RA (radians)
      }
      \sstsubsection{
         dh = double (Given)
      }{
         Hipparcos Dec (radians)
      }
      \sstsubsection{
         epoch = double (Given)
      }{
         Julian epoch (TDB)
      }
      \sstsubsection{
         r5 = double $*$ (Returned)
      }{
         RA (radians, FK5, equinox J2000, epoch \texttt{"} epoch\texttt{"} )
      }
      \sstsubsection{
         d5 = double $*$ (Returned)
      }{
         Dec (radians, FK5, equinox J2000, epoch \texttt{"} epoch\texttt{"} )
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraEpj2jd and eraHfk5z(). See SOFA/ERFA documentation for details.
      }
   }
}
\sstroutine{
   palRefcoq
}{
   Determine the constants A and B in the atmospheric refraction model
}{
   \sstdescription{
      Determine the constants A and B in the atmospheric refraction
      model dZ = A tan Z $+$ B tan$*$$*$3 Z.  This is a fast alternative
      to the palRefco routine.

      Z is the \texttt{"} observed\texttt{"}  zenith distance (i.e. affected by refraction)
      and dZ is what to add to Z to give the \texttt{"} topocentric\texttt{"}  (i.e. in vacuo)
      zenith distance.
   }
   \sstinvocation{
      palRefcoq( double tdk, double pmb, double rh, double wl,
                 double $*$refa, double $*$refb );
   }
   \sstarguments{
      \sstsubsection{
         tdk = double (Given)
      }{
         Ambient temperature at the observer (K)
      }
      \sstsubsection{
         pmb = double (Given)
      }{
         Pressure at the observer (millibar)
      }
      \sstsubsection{
         rh =  double (Given)
      }{
         Relative humidity at the observer (range 0-1)
      }
      \sstsubsection{
         wl =  double (Given)
      }{
         Effective wavelength of the source (micrometre).
         Radio refraction is chosen by specifying wl $>$ 100 micrometres.
      }
      \sstsubsection{
         refa = double $*$ (Returned)
      }{
         tan Z coefficient (radian)
      }
      \sstsubsection{
         refb = double $*$ (Returned)
      }{
         tan$*$$*$3 Z coefficient (radian)
      }
   }
   \sstnotes{
      \sstitemlist{

         \sstitem
          Uses eraRefco(). See SOFA/ERFA documentation for details.

         \sstitem
          Note that the SOFA/ERFA routine uses different order of
           of arguments and uses deg C rather than K.
      }
   }
}


% ? End of main text
\end{document}