source: trunk/MagicSoft/slalib/el2ue.c

#include "slalib.h"
#include "slamac.h"
void slaEl2ue ( double date, int jform, double epoch, double orbinc,
                double anode, double perih, double aorq, double e,
                double aorl, double dm, double u[], int *jstat )
/*
**  - - - - - - - - -
**   s l a E l 2 u e
**  - - - - - - - - -
**
**  Transform conventional osculating orbital elements into "universal" form.
**
**  Given:
**     date    double     epoch (TT MJD) of osculation (Note 3)
**     jform   int        choice of element set (1-3, Note 6)
**     epoch   double     epoch (TT MJD) of the elements
**     orbinc  double     inclination (radians)
**     anode   double     longitude of the ascending node (radians)
**     perih   double     longitude or argument of perihelion (radians)
**     aorq    double     mean distance or perihelion distance (AU)
**     e       double     eccentricity
**     aorl    double     mean anomaly or longitude (radians, jform=1,2 only)
**     dm      double     daily motion (radians, jform=1 only)
**
**  Returned:
**     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
**
**     jstat   int*       status:  0 = OK
**                                -1 = illegal jform
**                                -2 = illegal e
**                                -3 = illegal aorq
**                                -4 = illegal dm
**                                -5 = numerical error
**
**  Called:  slaUe2pv, slaPv2ue
**
**  Notes
**
**  1  The "universal" 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 "universal eccentric anomaly" at a given
**     date and (v) that date.
**
**  2  The companion routine is slaUe2pv.  This takes the set of numbers
**     that the present routine outputs and uses them to derive the
**     object's position and velocity.  A single prediction requires one
**     call to the present routine followed by one call to slaUe2pv;
**     for convenience, the two calls are packaged as the routine
**     slaPlanel.  Multiple predictions may be made by again calling the
**     present routine once, but then calling slaUe2pv multiple times,
**     which is faster than multiple calls to slaPlanel.
**
**  3  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).
**
**  4  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.
**
**  5  The universal elements returned in the array u are in canonical
**     units (solar masses, AU and canonical days).
**
**  6  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)
**
**  7  Unused elements (dm for jform=2, aorl and dm for jform=3) are
**     not accessed.
**
**  8  The algorithm was originally adapted from the EPHSLA program of
**     D.H.P.Jones (private communication, 1996).  The method is based on
**     Stumpff's Universal Variables.
**
**  Reference:  Everhart, E. & Pitkin, E.T., Am.J.Phys. 51, 712, 1983.
**
**  Last revision:   18 March 1999
**
**  Copyright P.T.Wallace.  All rights reserved.
*/

/* Gaussian gravitational constant (exact) */
#define GCON 0.01720209895

/* Sin and cos of J2000 mean obliquity (IAU 1976) */
#define SE 0.3977771559319137
#define CE 0.9174820620691818

{
   int j;
   double pht, argph, q, w, cm, alpha, phs, sw, cw, si, ci, so, co,
          x, y, z, px, py, pz, vx, vy, vz, dt, fc, fp, psi, ul[13], pv[6];



/* Validate arguments. */
   if ( jform < 1 || jform > 3 ) {
      *jstat = -1;
      return;
   }
   if ( e < 0.0 || e > 10.0 || ( e >= 1.0 && jform != 3 ) ) {
      *jstat = -2;
      return;
   }
   if ( aorq <= 0.0 ) {
      *jstat = -3;
      return;
   }
   if ( jform == 1 && dm <= 0.0 ) {
      *jstat = -4;
      return;
   }

/*
** Transform elements into standard form:
**
** pht   = epoch of perihelion passage
** argph = argument of perihelion (little omega)
** q     = perihelion distance (q)
** cm    = combined mass, M+m (mu)
*/

   if ( jform == 1 ) {
      pht = epoch - ( aorl - perih ) / dm;
      argph = perih - anode;
      q = aorq * ( 1.0 - e );
      w = dm / GCON;
      cm =  w * w * aorq * aorq * aorq;
   } else if ( jform == 2 ) {
      pht = epoch - aorl * sqrt ( aorq * aorq * aorq ) / GCON;
      argph = perih;
      q = aorq * ( 1.0 - e );
      cm = 1.0;
   } else if ( jform == 3 ) {
      pht = epoch;
      argph = perih;
      q = aorq;
      cm = 1.0;
   }

/*
** The universal variable alpha.  This is proportional to the total
** energy of the orbit:  -ve for an ellipse, zero for a parabola,
** +ve for a hyperbola.
*/

   alpha = cm * ( e - 1.0 ) / q;

/* Speed at perihelion. */

   phs = sqrt ( alpha + 2.0 * cm / q );

/*
** In a Cartesian coordinate system which has the x-axis pointing
** to perihelion and the z-axis normal to the orbit (such that the
** object orbits counter-clockwise as seen from +ve z), the
** perihelion position and velocity vectors are:
**
**   position   [Q,0,0]
**   velocity   [0,phs,0]
**
** To express the results in J2000 equatorial coordinates we make a
** series of four rotations of the Cartesian axes:
**
**           axis      Euler angle
**
**     1      z        argument of perihelion (little omega)
**     2      x        inclination (i)
**     3      z        longitude of the ascending node (big omega)
**     4      x        J2000 obliquity (epsilon)
**
** In each case the rotation is clockwise as seen from the +ve end
** of the axis concerned.
*/

/* Functions of the Euler angles. */
   sw = sin ( argph );
   cw = cos ( argph );
   si = sin ( orbinc );
   ci = cos ( orbinc );
   so = sin ( anode );
   co = cos ( anode );

/* Position at perihelion (AU). */
   x = q * cw;
   y = q * sw;
   z = y * si;
   y = y * ci;
   px = x * co - y * so;
   y = x * so + y * co;
   py = y * CE - z * SE;
   pz = y * SE + z * CE;

/* Velocity at perihelion (AU per canonical day). */
   x = - phs * sw;
   y = phs * cw;
   z = y * si;
   y = y * ci;
   vx = x * co - y * so;
   y = x * so + y * co;
   vy = y * CE - z * SE;
   vz = y * SE + z * CE;

/* Time from perihelion to date (in Canonical Days: a canonical */
/* day is 58.1324409... days, defined as 1/GCON).               */

   dt = ( date - pht ) * GCON;

/* First Approximation to the Universal Eccentric Anomaly, psi, */
/* based on the circle (fc) and parabola (fp) values.           */
   fc = dt / q;
   w = pow ( 3.0 * dt + sqrt ( 9.0 * dt * dt + 8.0 * q * q * q ),
             1.0 / 3.0 );
   fp = w - 2.0 * q / w;
   psi = ( 1.0 - e ) * fc + e * fp;

/* Assemble local copy of element set. */
   ul[0] = cm;
   ul[1] = alpha;
   ul[2] = pht;
   ul[3] = px;
   ul[4] = py;
   ul[5] = pz;
   ul[6] = vx;
   ul[7] = vy;
   ul[8] = vz;
   ul[9] = q;
   ul[10] = 0.0;
   ul[11] = date;
   ul[12] = psi;

/* Predict position+velocity at epoch of osculation. */
   slaUe2pv ( date, ul, pv, &j );
   if ( j ) {
      *jstat = -5;
      return;
   }

/* Convert back to universal elements. */
   slaPv2ue ( pv, date, cm - 1.0, u, &j );
   if ( j ) {
      *jstat = -5;
      return;
   }

/* OK exit. */
   *jstat = 0;

}