#include "slalib.h"
#include "slamac.h"
void slaPlantu ( double date, double elong, double phi, double u[],
                 double *ra, double *dec, double *r, int *jstat)
/*
**  - - - - - - - - - -
**   s l a P l a n t u
**  - - - - - - - - - -
**
**  Topocentric apparent RA,Dec of a Solar-System object whose
**  heliocentric universal elements are known.
**
**  Given:
**     date    double     TT MJD of observation (JD - 2400000.5)
**     elong   double     observer's east longitude (radians)
**     phi     double     observer's geodetic latitude (radians)
**     u       double[13] 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
**
**  Returned:
**     ra,dec  double*    RA, Dec (topocentric apparent, radians)
**     r       double*    distance from observer (AU)
**     jstat   int*       status:  0 = OK
**                                -1 = radius vector zero
**                                -2 = failed to converge
**
**  Notes:
**
**  1  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).
**
**  2  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 slaEvp and slaUe2pv.
**
**  3  The "universal" 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 "universal eccentric anomaly" at a given
**     date and (v) that date.
**
**  4  The universal elements are with respect to the J2000 equator and
**     equinox.
**
**  Called: slaEvp, slaUe2pv, slaGmst, slaDt, slaEpj, slaPrenut,
**          slaDmxv, slaPvobs, slaDcc2s, slaDranrm
**
**     1  Sterne, Theodore E., "An Introduction to Celestial Mechanics",
**        Interscience Publishers Inc., 1960.  Section 6.7, p199.
**
**     2  Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
**
**  Last revision:   27 December 2002
**
**  Copyright P.T.Wallace.  All rights reserved.
*/

/* Light time for unit distance (sec) */
#define TAU 499.004782

{
   int i;
   double dvb[3], dpb[3], vsg[6], vsp[6], v[6], rmat[3][3], vgp[6], stl,
          vgo[6], dx, dy, dz, d, tl;



/* Sun to geocentre (J2000). */
   slaEvp ( date, 2000.0, dvb, dpb, &vsg[3], &vsg[0] );

/* Sun to planet (J2000). */
   slaUe2pv ( date, u, vsp, jstat );

/* Geocentre to planet (J2000). */
   for ( i = 0; i < 6; i++ ) {
      v[i] = vsp[i] - vsg[i];
   }

/* Precession and nutation to date. */
   slaPrenut ( 2000.0, date, rmat );
   slaDmxv ( rmat, v, vgp );
   slaDmxv ( rmat, &v[3], &vgp[3] );

/* Geocentre to observer (date). */
   stl = slaGmst ( date - slaDt ( slaEpj ( date ) ) / 86400.0 ) + elong;
   slaPvobs ( phi, 0.0, stl, vgo );

/* Observer to planet (date). */
   for ( i = 0; i < 6; i++ ) {
      v[i] = vgp[i] - vgo[i];
   }

/* Geometric distance (AU). */
   dx = v[0];
   dy = v[1];
   dz = v[2];
   d = sqrt ( dx * dx + dy * dy + dz * dz );

/* Light time (sec). */
   tl = TAU * d;

/* Correct position for planetary aberration. */
   for ( i = 0; i < 3; i++ ) {
      v[i] = v[i] - tl * v[i+3];
   }

/* To RA,Dec. */
   slaDcc2s ( v, ra, dec );
   *ra = slaDranrm ( *ra );
   *r = d;

}
