Table of contents
CSPICE_EQNCPV computes 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).
Given:
et the scalar double precision ephemeris time, expressed as seconds
past J2000 TDB, at which the state of the target body is to be
computed.
help, et
DOUBLE = Scalar
`et' has units of TDB seconds.
epoch the scalar double precision epoch of the equinoctial elements in
seconds past the J2000 epoch.
help, epoch
DOUBLE = Scalar
`epoch' has units of TDB seconds.
eqel a double precision 9-vector containing 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.
help, eqel
DOUBLE = Array[9]
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[0] is the semi-major axis (A) of the orbit in
km.
eqel[1] is the value of H at the specified epoch.
( E*sin(argp+node) ).
eqel[2] is the value of K at the specified epoch
( E*cos(argp+node) ).
eqel[3] is the mean longitude (mean0+argp+node) at
the epoch of the elements measured in
radians.
eqel[4] is the value of P (tan(inc/2)*sin(node)) at
the specified epoch.
eqel[5] is the value of Q (tan(inc/2)*cos(node)) at
the specified epoch.
eqel[6] 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[7] 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[8] 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 the scalar double precision Right Ascension of the pole of the
reference plane with respect to some inertial frame.
help, rapol
DOUBLE = Scalar
`rapol' has units of radians.
decpol the scalar double precision Declination of the pole of the
reference plane with respect to some inertial frame.
help, decpol
DOUBLE = Scalar
`decpol' has units of radians.
the call:
cspice_eqncpv, et, epoch, eqel, rapol, decpol, state
returns:
state the double precision 6-vector describing the state of the object
defined by `eqel' relative to the inertial frame used to define
`rapol' and `decpol'.
help, state
DOUBLE = Array[6]
Units are in km and km/sec.
None.
Any numerical results shown for this example may differ between
platforms as the results depend on the SPICE kernels used as input
and the machine specific arithmetic implementation.
1) Compute a state vector from a set of equinoctial elements.
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 computes the
state of the object at an arbitrary epoch. The code below
illustrates how to 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.
domega Rate of change of argument of periapse in
radians/second.
epoch is the epoch of the elements in seconds past
the J2000 epoch.
Example code begins here.
PRO eqncpv_ex1
;;
;; eqel[0] = a
;; eqel[1] = ecc * sin( omega + node )
;; eqel[2] = ecc * cos( omega + node )
;;
;; eqel[3] = m + omega + node
;;
;; eqel[4] = tan(inc/2.0) * sin(node)
;; eqel[5] = tan(inc/2.0) * cos(node)
;;
;; eqel[6] = domega
;; eqel[7] = domega + dmdt + dnode
;; eqel[8] = dnode
;;
;; In this case, the rates of node and argument of
;; periapse are zero and the pole of the central
;; frame is aligned with the pole of an inertial frame.
;;
p = 1.0d4
gm = 398600.436d
ecc = 0.1d
a = p/( 1.d - ecc )
dmdt = sqrt ( gm / a ) / a
omega = 30.d * cspice_rpd()
node = 15.d * cspice_rpd()
inc = 10.d * cspice_rpd()
m = 45.d * cspice_rpd()
epoch = -100000000.d
dnode = 0.d
domega = 0.d
eqel = [ a, $
ecc*sin(omega+node), $
ecc*cos(omega+node), $
m + omega + node, $
tan(inc/2.d)*sin(node), $
tan(inc/2.d)*cos(node), $
domega, $
domega + dmdt + dnode, $
dnode ]
rapol = -cspice_halfpi()
decpol = cspice_halfpi()
et = epoch - 10000.d0
for i = 1, 10 do begin
et = et + 250.d
print
cspice_eqncpv, et, epoch, eqel, rapol, decpol, state
print, 'Pos = ',state[0:2]
print, 'Vel = ',state[3:5]
endfor
END
When this program was executed on a Mac/Intel/IDL8.x/64-bit
platform, the output was:
Pos = -10732.167 3902.5058 1154.4516
Vel = -2.5407669 -5.1522692 -0.76157581
Pos = -11278.383 2586.1799 955.18410
Vel = -1.8271564 -5.3629158 -0.83001977
Pos = -11645.295 1228.6124 740.70957
Vel = -1.1080964 -5.4828109 -0.88325573
Pos = -11832.800 -147.99098 514.80525
Vel = -0.39342066 -5.5159047 -0.92150772
Pos = -11843.089 -1522.4698 281.17526
Vel = 0.30828845 -5.4665647 -0.94512794
Pos = -11680.365 -2874.7848 43.424394
Vel = 0.98951988 -5.3393639 -0.95455246
Pos = -11350.590 -4186.0498 -194.95853
Vel = 1.6436489 -5.1389380 -0.95026851
Pos = -10861.293 -5438.5362 -430.61041
Vel = 2.2647587 -4.8698990 -0.93279157
Pos = -10221.411 -6615.6606 -660.29899
Vel = 2.8474759 -4.5367943 -0.90265092
Pos = -9441.1703 -7701.9679 -880.92519
Vel = 3.3868216 -4.1441030 -0.86038222
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.
1) If the eccentricity corresponding to the input elements is
greater than 0.9, the error SPICE(ECCOUTOFRANGE) is signaled
by a routine in the call tree of this routine.
2) If the semi-major axis of the elements is non-positive, the
error SPICE(BADSEMIAXIS) is signaled by a routine in the call
tree of this routine.
3) If any of the input arguments, `et', `epoch', `eqel', `rapol'
or `decpol', is undefined, an error is signaled by the IDL
error handling system.
4) If any of the input arguments, `et', `epoch', `eqel', `rapol'
or `decpol', is not of the expected type, or it does not have
the expected dimensions and size, an error is signaled by the
Icy interface.
5) If the output argument `state' is not a named variable, an
error is signaled by the Icy interface.
None.
1) The equinoctial elements used by this routine are taken
from "Tangent" formulation of equinoctial elements
P = tan(inclination/1) * sin(R.A. of ascending node)
Q = tan(inclination/1) * 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 / sqrt ( 1.0 - P*P - Q*Q )
Q = Q / sqrt ( 1.0 - P*P - Q*Q )
This will convert the Sine formulation to the Tangent
formulation.
ICY.REQ
SPK.REQ
[1] W. Owen and R. Vaughan, "Optical Navigation Program
Mathematical Models," JPL Engineering Memorandum 314-513,
August 9, 1991.
J. Diaz del Rio (ODC Space)
E.D. Wright (JPL)
-Icy Version 1.0.1, 31-MAY-2021 (JDR)
Edited the header to comply with NAIF standard. Added
example's problem statement.
Added -Parameters, -Exceptions, -Files, -Restrictions,
-Literature_References and -Author_and_Institution sections.
Removed reference to the routine's corresponding CSPICE header from
-Abstract section.
Added arguments' type and size information in the -I/O section.
-Icy Version 1.0.0, 22-NOV-2011 (EDW)
Compute a state from equinoctial elements
|