Procedure

     EQNCPV (Equinoctial Elements to position and velocity)

     SUBROUTINE EQNCPV ( ET, EPOCH, EQEL, RAPOL, DECPOL, STATE )

Abstract

     Compute the state (position and velocity) of an object whose
     trajectory is described via equinoctial elements relative to some
     fixed plane (usually the equatorial plane of some planet).

Required_Reading

     SPK

Keywords

     EPHEMERIS

Declarations

     IMPLICIT NONE

     DOUBLE PRECISION      ET
     DOUBLE PRECISION      EPOCH
     DOUBLE PRECISION      EQEL ( 9 )
     DOUBLE PRECISION      RAPOL
     DOUBLE PRECISION      DECPOL
     DOUBLE PRECISION      STATE ( 6 )

Brief_I/O

     VARIABLE  I/O  DESCRIPTION
     --------  ---  --------------------------------------------------
     ET         I   Epoch in seconds past J2000 to find state
     EPOCH      I   Epoch of elements in seconds past J2000
     EQEL       I   Array of equinoctial elements
     RAPOL      I   Right Ascension of the pole of the reference plane
     DECPOL     I   Declination of the pole of the reference plane
     STATE      O   State of the object described by EQEL.

Detailed_Input

     ET       is the epoch (ephemeris time) at which the state
              of the target body is to be computed. ET is measured
              in seconds past the J2000 epoch.

     EPOCH    is the epoch of the equinoctial elements in seconds
              past the J2000 epoch.

     EQEL     is an array of 9 double precision numbers that are the
              equinoctial elements for some orbit expressed relative to
              the equatorial frame of the central body defined as

              -  The Z-axis of the equatorial frame is the direction
                 of the pole of the central body relative to some
                 inertial frame;

              -  The X-axis is given by the cross product of the Z-axis
                 of the inertial frame with the direction of the pole
                 of the central body; and

              -  The Y-axis completes a right handed frame.

              If the X-axis of the equatorial frame is aligned with the
              X-axis of the inertial frame, then the X-axis of the
              equatorial frame will be located at 90 degrees + RAPOL in
              the inertial frame.

              The specific arrangement of the elements is spelled out
              below:

                 EQEL(1)   is the semi-major axis (A) of the orbit in
                           km.

                 EQEL(2)   is the value of H at the specified epoch.
                           ( E*SIN(ARGP+NODE) ).

                 EQEL(3)   is the value of K at the specified epoch
                           ( E*COS(ARGP+NODE) ).

                 EQEL(4)   is the mean longitude (MEAN0+ARGP+NODE) at
                           the epoch of the elements measured in
                           radians.

                 EQEL(5)   is the value of P (TAN(INC/2)*SIN(NODE)) at
                           the specified epoch.

                 EQEL(6)   is the value of Q (TAN(INC/2)*COS(NODE)) at
                           the specified epoch.

                 EQEL(7)   is the rate of the longitude of periapse
                           (dARGP/dt + dNODE/dt ) at the epoch of
                           the elements. This rate is assumed to hold
                           for all time. The rate is measured in
                           radians per second.

                 EQEL(8)   is the derivative of the mean longitude
                           ( dM/dt + dARGP/dt + dNODE/dt ). This
                           rate is assumed to be constant and is
                           measured in radians/second.

                 EQEL(9)   is the rate of the longitude of the
                           ascending node ( dNODE/dt). This rate is
                           measured in radians per second.

              where

                 INC       is the inclination of the orbit,

                 ARGP      is the argument of periapse,

                 NODE      is longitude of the ascending node, and

                 E         is eccentricity of the orbit.

     RAPOL    is the Right Ascension of the pole of the reference plane
              with respect to some inertial frame (measured in
              radians).

     DECPOL   is the Declination of the pole of the reference plane
              with respect to some inertial frame (measured in
              radians).

Detailed_Output

     STATE    is the state of the object described by EQEL relative to
              the inertial frame used to define RAPOL and DECPOL. Units
              are in km and km/sec.

Parameters

     None.

Exceptions

     1)  If the eccentricity corresponding to the input elements is
         greater than 0.9, the error SPICE(ECCOUTOFRANGE) is signaled.

     2)  If the semi-major axis of the elements is non-positive, the
         error SPICE(BADSEMIAXIS) is signaled.

Files

     None.

Particulars

     This routine evaluates the input equinoctial elements for
     the specified epoch and return the corresponding state.

     This routine was adapted from a routine provided by
     Bob Jacobson of the Planetary Dynamics Group of
     the Navigation and Flight Mechanics Section at JPL.

Examples

     Suppose you have classical elements and rates of
     change of the ascending node and argument of periapse
     for some satellite of the earth.

     By transforming the classical elements
     this routine can be used to compute the state of the
     object at an arbitrary epoch. The code below illustrates
     how you might do this.

     The table below illustrates the meanings of the various
     variables used in the discussion below.

           Variable     Meaning
           --------     ----------------------------------
           A            Semi-major axis in km
           ECC          Eccentricity of orbit
           INC          Inclination of orbit
           NODE         Longitude of the ascending node at epoch
           OMEGA        Argument of periapse at epoch
           M            Mean anomaly at epoch
           DMDT         Mean anomaly rate in radians/second
           DNODE        Rate of change of longitude of ascending node
                        in radians/second
           DARGP        Rate of change of argument of periapse in
                        radians/second
           EPOCH        is the epoch of the elements in seconds past
                        the J2000 epoch.


        EQEL(1) = A
        EQEL(2) = ECC * DSIN ( OMEGA + NODE )
        EQEL(3) = ECC * DCOS ( OMEGA + NODE )

        EQEL(4) = M + OMEGA + NODE

        EQEL(5) = TAN(INC/2.0D0) * DSIN(NODE)
        EQEL(6) = TAN(INC/2.0D0) * DCOS(NODE)

        EQEL(7) = DARGP
        EQEL(8) = DARGP + DMDT + DNODE
        EQEL(9) = DNODE


        We shall compute the state of the satellite in the
        pole and equator reference system.

        RAPOL   = -HALFPI()
        DECPOL  =  HALFPI()


        Now compute the state at the desired epoch ET.

        CALL EQNCPV ( ET, EPOCH, EQEL, RAPOL, DECPOL, STATE )

Restrictions

     1)  The equinoctial elements used by this routine are taken
         from  "Tangent" formulation of equinoctial elements

            p = tan(inclination/2) * sin(R.A. of ascending node)
            q = tan(inclination/2) * cos(R.A. of ascending node)

         Other formulations use Sine instead of Tangent. We shall
         call these the "Sine" formulations.

            p = sin(inclination/2) * sin(R.A. of ascending node)
            q = sin(inclination/2) * cos(R.A. of ascending node)

         If you have equinoctial elements from this alternative
         formulation you should replace p and q  by the
         expressions below.

            P = P / DSQRT ( 1.0D0 - P*P - Q*Q )
            Q = Q / DSQRT ( 1.0D0 - P*P - Q*Q )

         This will convert the Sine formulation to the Tangent
         formulation.

Literature_References

     [1]  W. Owen and R. Vaughan, "Optical Navigation Program
          Mathematical Models," JPL Engineering Memorandum 314-513,
          August 9, 1991.

Author_and_Institution

     J. Diaz del Rio    (ODC Space)
     R.A. Jacobson      (JPL)
     B.V. Semenov       (JPL)
     W.L. Taber         (JPL)

Version

    SPICELIB Version 1.0.3, 14-APR-2021 (JDR)

        Edited the header to comply with NAIF standard. Added SPK
        required reading.

    SPICELIB Version 1.0.2, 18-MAY-2010 (BVS)

        Removed "C$" marker from text in the header.

    SPICELIB Version 1.0.1, 31-JAN-2008 (BVS)

        Removed non-standard header section heading
        'Declarations_of_external_functions'.

    SPICELIB Version 1.0.0, 08-JAN-1997 (WLT) (RAJ)