function aberat (pos1, ve, tlight, pos2)

` this function corrects position vector for aberration of light.
` algorithm includes relativistic terms.  see murray (1981)
` mon. notices royal ast. society 195, 639-648.

` input

`  pos1   = position vector, referred to origin at center of
`           mass of the earth, components in au
`  ve     = velocity vector of center of mass of the earth,
`           referred to origin at solar system barycenter,
`           components in au/day
`  tlight = light time from body to earth in days
`           if tlight = 0, this function will compute tlight

` output

`   pos2 = position vector, referred to origin at center of
`          mass of the earth, corrected for aberration,
`          components in au

--------------------------------------------------

function apparent(tjd, ujd, l, n, topo, glon, glat, ht, ra, dec, dis)

` this function computes the apparent geocentric or topocentric place
` of a planet or other solar system body. rectangular coordinates of
` solar system bodies are obtained from an slp96 binary data file.

` input

`  tjd  = tdt julian date for apparent geocentric place
`  ujd  = ut1 julian date for apparent topocentric place
`  l    = body identification number for desired planet
`  n    = body identification number for the earth
`  topo = type of apparent place calculation
`       = 0 => geocentric
`       = 1 => topocentric
`  glon = geodetic longitude of observer (east +, degrees)
`  glat = geodetic latitude of observer (north +, degrees)
`  ht   = height of observer (meters)

` output

`  ra  = apparent geocentric or topocentric right ascension,
`        referred to true equator and equinox of date (hours)
`  dec = apparent geocentric or topocentric declination,
`        referred to true equator and equinox of date (degrees)
`  dis = true distance from earth to planet (astronomical units)

--------------------------------------------------

function apstar (tjd, ujd, topo, glon, glat, ht, ram, decm, pmra, pmdec, \
                 parlax, radvel, ra, dec)

` this function computes the geocentric or topocentric apparent place
` of a star, given its mean place, proper motion, parallax, and radial
` velocity for j2000.0.

` input

`  tjd    = tdt julian date for apparent place
`  ujd    = ut1 julian date for apparent topocentric place
`  topo   = type of apparent place calculation
`         = 0 ==> geocentric
`         = 1 ==> topocentric
`  glon   = geodetic longitude of observer (east +, degrees)
`  glat   = geodetic latitude of observer (north +, degrees)
`  ht     = height of observer (meters)
`  ram    = mean right ascension j2000.0 (hours)
`  decm   = mean declination j2000.0 (degrees)
`  pmra   = proper motion in ra (seconds of time/julian century)
`  pmdec  = proper motion in dec (seconds of arc/julian century)
`  parlax = parallax (seconds of arc)
`  radvel = radial velocity (kilometers per second)

` output

`  ra  = apparent geocentric or topocentric right ascension,
`        referred to true equator and equinox of date (hours)
`  dec = apparent geocentric or topocentric declination,
`        referred to true equator and equinox of date (degrees)

--------------------------------------------------

function earth (jdate, rearth)

` true-of-date heliocentric, ecliptic
` position vector of the Earth

` input

`  jdate = julian day

` output

`  rearth = position vector of the earth (km)

--------------------------------------------------

function elp2000(tjd, oev)

` osculating orbital elements of the moon

` mean ecliptic and mean equinox of date

` input

`  tjd = tdt julian date

` output

`  oev[1]  = semimajor axis (kilometers)
`  oev[2]  = orbital eccentricity (non-dimensional)
`            (0 <= eccentricity < 1)
`  oev[3]  = orbital inclination (radians)
`            (0 <= inclination <= pi)
`  oev[4]  = argument of perigee (radians)
`            (0 <= argument of perigee <= 2 pi)
`  oev[5]  = right ascension of ascending node (radians)
`            (0 <= raan <= 2 pi)
`  oev[6]  = true anomaly (radians)
`            (0 <= true anomaly <= 2 pi)

--------------------------------------------------

function etilt (tjd, oblm, oblt, eqeq, dpsi, deps)

` this function computes quantities related to the orientation
` of the earth's rotation axis at julian date tjd

` input

`  tjd = tdb julian date for orientation parameters

` output

`  oblm = mean obliquity of the ecliptic at date tjd (degrees)
`  oblt = true obliquity of the ecliptic at date tjd (degrees)
`  eqeq = equation of the equinoxes at date tjd (arc seconds)
`  dpsi = nutation in longitude at date tjd (arc seconds)
`  deps = nutation in obliquity at date tjd (arc seconds)

--------------------------------------------------

function gast1 (tjdh, tjdl, k, gst)

` this function computes the greenwich sidereal time
` (either mean or apparent) at julian date tjdh + tjdl

` input

`  tjdh = julian date, high-order part

`  tjdl = julian date, low-order part

`         julian date may be split at any point, but for
`         highest precision, set tjdh to be the integral part of
`         the julian date, and set tjdl to be the fractional part

`  k    = time selection code
`         set k=0 for greenwich mean sidereal time
`         set k=1 for greenwich apparent sidereal time

` output

`  gst = greenwich (mean or apparent) sidereal time in hours

--------------------------------------------------

function gast2 (jdate, gst)

` Greenwich apparent sidereal time

` input

`  jdate = Julian date

` output

`  gst = Greenwich apparent sidereal time (radians)
`        (0 <= gst <= 2 pi)

--------------------------------------------------

function geocen (pos1, pe, pos2, tlight)

` this function moves the origin of coordinates from the
` barycenter of the solar system to the center of mass of the
` earth, i.e., this function corrects for parallax.

` input

`  pos1 = position vector, referred to origin at solar system
`         barycenter (au)
`  pe   = position vector of center of mass of the earth,
`         referred to origin at solar system barycenter (au)

` output

`   pos2   = position vector, referred to origin at center of
`            mass of the earth (au)
`   tlight = light time from body to earth in days

--------------------------------------------------

function jupiter (jdate, rjupiter)

` true-of-date heliocentric, ecliptic
` position vector of Jupiter

` input

`  jdate = julian day

` output

`  rjupiter = position vector of Jupiter (km)

--------------------------------------------------

function mars (jdate, rmars)

` true-of-date heliocentric, ecliptic
` position vector of Mars

` input

`  jdate = julian day

` output

`  rmars = position vector of mars (km)

--------------------------------------------------

function moon (jed, rmoon)

` geocentric position of the moon

` input

`  jed = julian ephemeris date

` output

`  rmoon = geocentric, equatorial position
`          vector of the moon (kilometers)

--------------------------------------------------

function nutate (tjd, pos1, pos2)

` this function nutates equatorial rectangular coordinates from
` mean equator and equinox of epoch to true equator and equinox of
` epoch. see pages 41-45 of the explanatory supplement to the ae.

` input

`  tjd  = tdb julian date of epoch
`  pos1 = position vector, geocentric equatorial rectangular
`         coordinates, referred to mean equator and equinox
`         of epoch

` output

`   pos2 = position vector, geocentric equatorial rectangular
`          coordinates, referred to true equator and equinox
`          of epoch (out)

` note:
` if tjd is negative, inverse nutation (true to mean) is applied

--------------------------------------------------

function nutation (jdate, dpsi, deps)

` nutation in longitude and obliquity

` IAU 1980 version

` input

`  jdate = julian date (tdb)

` output

`  dpsi = nutation in longitude (arc seconds)
`  deps = nutation in obliquity (arc seconds)

--------------------------------------------------

function pluto(jdate, rpluto)

` heliocentric coordinates of pluto

` heliocentric position vector of Pluto
` ecliptic and equinox of J2000

` input

`  jdate = julian ephemeris date

` output

`  rpluto = position vector of Pluto (km)

--------------------------------------------------

function precess (tjd1, pos1, tjd2, pos2)

` this function precesses equatorial rectangular coordinates from
` one epoch to another.  the coordinates are referred to the mean
` equator and equinox of the two respective epochs. see pages 30-34
` of the explanatory supplement to the ae, lieske, et al. (1977)
` astronomy and astrophysics 58, 1-16, and lieske (1979) astronomy
` and astrophysics 73, 282-284.

` input

`  tjd1 = tdb julian date of first epoch

`  pos1 = position vector, geocentric equatorial rectangular
`         coordinates, referred to mean equator and equinox of
`         first epoch
`  tjd2 = tdb julian date of second epoch

` output

`  pos2 = position vector, geocentric equatorial rectangular
`         coordinates, referred to mean equator and equinox of
`         second epoch

--------------------------------------------------

function propmo (tjd1, pos1, vel1, tjd2, pos2)

` this function applies proper motion, including
` foreshortening effects, to a star's position.

` input

`  tjd1 = tdb julian date of first epoch
`  pos1 = position vector at first epoch
`  vel1 = velocity vector at first epoch
`  tjd2 = tdb julian date of second epoch

` output

`  pos2 = position vector at second epoch

--------------------------------------------------

function readstar(fname, sname, ram, decm, pmra, pmdec, parlax, radvel)

` read star data file

` input

`  fname = name of star data file

` output

`  sname  = star name
`  ram    = J2000 right ascension (hours)
`  decm   = J2000 declination (degrees)
`  pmra   = J2000 proper motion in right ascension
`           (seconds/Julian century)
`  pmdec  = J2000 proper motion in declination
`           (arcseconds/Julian century)
`  parlax = parallax (arcseconds)
`  radvel = radial velocity (kilometers/second)

--------------------------------------------------

function saturn (jdate, rsaturn)

` true-of-date heliocentric, ecliptic
` position vector of Saturn

` input

`  jdate = julian day

` output

`  rsaturn = position vector of Saturn (km)

--------------------------------------------------

function slp96(tjd, ibody, icent, ipv, iframe, icoord, result, ierr)

` solar, lunar and planetary coordinates j2000

` input

`  tjd       julian date tdb

`  ibody     body index

`            ibody=01 : mercury       ibody=07 : uranus
`            ibody=02 : venus         ibody=08 : neptune
`            ibody=03 : e-m baryc.    ibody=09 : void
`            ibody=04 : mars          ibody=10 : moon
`            ibody=05 : jupiter       ibody=11 : sun
`            ibody=06 : saturn        ibody=12 : earth

`  icent     frame center index

`            icent=00 : barycenter of solar system
`            icent=01 : mercury       icent=07 : uranus
`            icent=02 : venus         icent=08 : neptune
`            icent=03 : e-m baryc.    icent=09 : void
`            icent=04 : mars          icent=10 : moon
`            icent=05 : jupiter       icent=11 : sun
`            icent=06 : saturn        icent=12 : earth

`  ipv       position-velocity index

`            ipv=1 : position.
`            ipv=2 : position and velocity.

`  iframe    frame index

`            iframe=1 : equinox and equator  j2000 (fk5)
`            iframe=2 : equinox and ecliptic j2000 (dynamical)

`  icoord    coordinates index

`            icoord=1 : rectangular coordinates
`            icoord=2 : spherical coordinates

`  fname     name of the slp96 data file

` output

`  result    results table

`            rectangular coordinates (icoord = 1)

`            position

`            result[1] : equatorial or ecliptic component x (au)
`            result[2] : equatorial or ecliptic component y (au)
`            result[3] : equatorial or ecliptic component z (au)

`            velocity

`            result[4] : equatorial or ecliptic component x (au/day)
`            result[5] : equatorial or ecliptic component y (au/day)
`            result[6] : equatorial or ecliptic component z (au/day)

`            spherical coordinates (icoord = 2)

`            position

`            result[1] : right ascension or longitude (radians)
`            result[2] : declination or latitude (radians)
`            result[3] : geometric distance (au)

`            velocity

`            result[4] : right ascension or longitude (rad/day)
`            result[5] : declination or latitude (rad/day)
`            result[6] : geometric distance (au/day)

`  ierr      error index

`            ierr=0  : no error
`            ierr=10 : file error
`            ierr=11 : date error (tjd)
`            ierr=12 : body error (ibody)
`            ierr=13 : frame center error (icent)
`            ierr=14 : position-velocity error (ipv)
`            ierr=15 : frame error (iframe)
`            ierr=16 : coordinates error (icoord)

--------------------------------------------------

function sun (jdate, rasc, decl, rsun)

` solar ephemeris

` input

`  jdate = julian day

` output

`  rasc = right ascension of the sun (radians)
`         (0 <= rasc <= 2 pi)
`  decl = declination of the sun (radians)
`         (-pi/2 <= decl <= pi/2)
`  rsun = eci position vector of the sun (km)

` note

`  coordinates are inertial, geocentric,
`  equatorial and true-of-date

--------------------------------------------------

function sun1 (jdate, rasc, decl, rsun)

` precision ephemeris of the Sun

` input

`  jdate = julian ephemeris date

` output

`  rasc = right ascension of the Sun (radians)
`         (0 <= rasc <= 2 pi)
`  decl = declination of the Sun (radians)
`         (-pi/2 <= decl <= pi/2)
`  rsun = eci position vector of the Sun (kilometers)

` note

`  coordinates are inertial, geocentric,
`  equatorial and true-of-date

--------------------------------------------------

function sunfld (pos1, pe, pos2)

` this function corrects position vector for the deflection
` of light in the gravitational field of the sun. see murray (1981)
` mon. notices royal ast. society 195, 639-648. this function valid
` for bodies within the solar system as well as for stars.

` input

`  pos1 = position vector, referred to origin at center of mass
`         of the earth, components in au
`  pe   = position vector of center of mass of the earth,
`         referred to origin at center of mass of the sun,
`         components in au

` output

`  pos2 = position vector, referred to origin at center of mass
`         of the earth, corrected for gravitational deflection,
`         components in au

--------------------------------------------------

function tdt2tdb(jdtdt, jdtdb)

` convert terrestrial dynamic time julian
` date to barycentric dynamic time julian date

` input

`  jdtdt = tdt julian date

` output

`  jdtdb = tdb julian date

--------------------------------------------------

function tdtimes (tdbjd, ttjd, secdif)

` this function computes the terrestrial time (tt)
` julian date corresponding to a barycentric dynamical
` time (tdb) julian date. expressions used in this
` version are approximations resulting in accuracies
` of about 20 microseconds. see explanatory supplement
` to the astronomical almanac, pp. 42-44 and 316.

` input

`  tdbjd = tdb julian date

` output

`  ttjd = tt julian date

`  secdif = difference tdbjd-ttjd, in seconds (out)

--------------------------------------------------

function terra (glon, glat, ht, st, posv, vel)

` this functions computes the position and velocity vectors of
` a terrestrial observer with respect to the center of the earth.

` input

`  glon = longitude of observer with respect to reference
`         meridian (east +) in degrees
`  glat = geodetic latitude (north +) of observer in degrees
`  ht   = height of observer in meters
`  st   = local apparent sidereal time at reference meridian
`         in hours

` output

`  posv = position vector of observer with respect to center
`         of earth, equatorial rectangular coordinates,
`         referred to true equator and equinox of date,
`         components in au
`  vel  = velocity vector of observer with respect to center
`         of earth, equatorial rectangular coordinates,
`         referred to true equator and equinox of date,
`         components in au/day

`  note: if reference meridian is greenwich and st=0, posv
`        is effectively referred to equator and greenwich

--------------------------------------------------

function utc2tdt(jdutc, jdtdt)

` convert utc julian date to tdt julian date

` input

`  jdutc = universal time julian date

` output

`  jdtdt = terrestrial dynamic time julian date

--------------------------------------------------

function venus (jdate, rvenus)

` true-of-date heliocentric, ecliptic
` position vector of venus

` input

`  jdate = julian day

` output

`  rvenus = position vector of venus (km)

--------------------------------------------------