Orbital elements¶
Ths computation of orbital element is either a clone of the NASA GMAT C++ code or of the Vallado algorithms (4th edition). The validation for the computation of all elements is listed at the bottom of this page, apart for the BPlane validation, which at the bottom of the BPlane section here.
API documentation available here for Orbit and here for OrbitDual.
Storage¶
Nyx stores all of the orbit information as a Cartesian state (units of kilometer and kilometer per second) because it is a nonsingular representation of an orbit. Furthermore, all propagation using OrbitalDynamics
, its acceleration models, and its force models is in Cartesian form for the same reason.
Initialization methods¶
An orbit may be initialized from its Cartesian position and velocity components (from 64 bit floating point values or from a vector of that type), from Keplerian orbital elements where the phasing parameter is the true anomaly, or from its geodesic elements (latitude, longitude and height compared to the reference ellipsoid of the frame^{1}). The initializer also includes an epoch (cf. Time MathSpec) and a frame (cf. Frame MathSpec).
As discussed above, the state is stored in Cartesian form. Hence, initializing in Keplerian form will trigger the conversion of the state from Keplerian orbital element into Cartesian orbital elements. This method does support hyperbolic, circular inclined, circular equatorial, and the common elliptical orbits. The eccentricity tolerance is set to \(1e^{11}\), i.e. if the eccentricity is below that number, then the orbit is considered circular and the appropriate conversions to Cartesian will be triggered. The algorithm implementation is available here, but to convince yourself that it works, probably best to check out the validation below.
Vallado's geodetic to Cartesian initializer is also implemented allowing the initialization of a state known only by its longitude, laltitude and height above the reference ellipsoid of that rocky celestial object.
Finally, one may also initialize a new orbit from a previous orbit via with_sma
, with_ecc
, with_inc
, with_raan
, with_aop
, with_ta
: these copy a previous orbit but change the requested Keplerian orbital element to the requested value. Analogous methods exist with add
instead of with
which will add the provided value to the current Keplerian element value. For example some_orbit.add_inc(10.0)
will create a new orbit where the inclination 10 degrees lower than some_orbit
's inclination.
Accessor methods¶
In the following, the position coordinates are referred to as \(x,y,z\) and the velocity components as \(v_x,v_y,v_z\). More generally, \(\gamma_x\), \(\gamma_y\), \(\gamma_z\) respectively refer to the x, y, and z component of the \(\gamma\) vector.
Argument of latitude (aol
)¶
The true longitude is an alternate phasing element. Nyx will return this angle between 0 and 360 degrees.

If the eccentricity is below the eccentricity tolerance:
\[ u = \lambda  \Omega \] 
Else:
\[ u = \omega + \nu \]
Argument of periapsis (aop
)¶
The argument of periapsis is computed as follows and is returned in degrees:

Rotate the momentum vector:
\[\begin{equation} \mathbf n = \begin{bmatrix} 0 \\ 0 \\ 1 \\ \end{bmatrix} \times \mathbf h \end{equation}\] 
Compute the AoP:
\[ \omega = \cos ^{1} \left(\frac{\mathbf n \cdot \mathbf e}{\mathbf n \times e}\right) \] 
Perform a quadrant check: if \(e_z\) is less then zero, return the AoP as \(2\pi  \omega\) instead.
Bplane (b_plane
)¶
Warning
This function will return an error (NyxError::NotHyperbolic
) if called on a nonhyperbolic orbit. In other words, to compute the Bplane of an orbit, convert it to the desired planetary inertial frame first.
The algorithm is identical to that outlined in the GMAT MathSpec, section 3.2.7, paraphrased here.
The \(\mathbf B\) vector is defined from the center of mass of the desired central body to the incoming hyperbolic asymptote: B is perpendicular to the incoming asymptote. We define \(\mathbf {\hat S}\) as a unit vector in the direction of the incoming asymptot and \(\mathbf{\hat T}\) as the unit vector perpendicular to \(\mathbf {\hat S}\) such that \(\mathbf {\hat T}\) lies in the \(xy\) plane of the frame in which is represented this orbit. \(\mathbf {\hat R}\) is a unit vector perpendicular to both \(\mathbf {\hat T}\) and \(\mathbf {\hat S}\).
We define \(B_T=\mathbf B \cdot \mathbf {\hat T}\) and \(B_R=\mathbf B \cdot \mathbf {\hat R}\). The method was adopted from work by Kizner

Compute the eccentricity unit vector, the momentum unit vector, and their cross product:
\[ \mathbf{\hat e} = \frac {\mathbf{e}}{e} \quad\quad \mathbf{\hat h} = \frac {\mathbf{h}}{h} \quad\quad \mathbf{\hat n} = \mathbf{\hat h}\times \mathbf{\hat e}\] 
Compute the \(\mathbf{\hat S}\), \(\mathbf B\), \(\mathbf{\hat T}\) and \(\mathbf{\hat R}\) vectors, where \(b\) is the semiminor axis.
\[ \mathbf{\hat S} = \frac {\mathbf{\hat e}} {e} + \sqrt{1  \left(\frac {1} {e}\right)^2} \mathbf{\hat n}\]\[ \mathbf{\hat T} = \frac {\mathbf{\hat S} \times \mathbf{\hat k}} {  \mathbf{\hat S} \times \mathbf{\hat k} } \]\[ \mathbf{B} = b \left(\sqrt{1 \left(\frac 1 e \right)^2 \mathbf{\hat e}}  \frac 1 e \mathbf{\hat n} \right)\]\[ \mathbf{\hat R} = \mathbf{\hat S} \times \mathbf{\hat T}\] 
Compute the B plane parameters
\[B_T=\mathbf B \cdot \mathbf {\hat T} \quad\quad B_R=\mathbf B \cdot \mathbf {\hat R}\]
Validation
The validation was done by following the Ex_LunarTransfer.script
GMAT example script, from which we've added the export of the B Plane parameters (B dot R, B dot T, B Vector Mag, B Vector Angle and C3).
The validation test is called val_b_plane_gmat
. It propagates the initial orbit of the script until Earth periapse and generates the trajectory.
Note: in that example, GMAT will use all of the force models available (high fidelity drag, SRP, all celestial objects) but we only use the point masses of Earth, Moon, Sun and Jupiter. This leads to a small error (less than 500 meters) for the B Plane. For each of the nine exported states (from some Earth apoapse to the Earth periapse), we check that the B plane and C3 parameters computed by Nyx are within 500 meters of error (due to force model differences) for positions and less than 0.1 millidegrees for B vector angle.
cargo test val_b_plane_gmat release  nocapture
B Vector Angle (angle
)¶
Returns the angle of the B plane, in degrees between \(]180;180]\).
B Vector magnitude (mag
)¶
Returns the magnitude of the B vector in kilometers.
B Plane frame from inertial frame (inertial_to_bplane
)¶
Returns the DCM to convert to this BPlane frame from the inertial frame (with identical central bodies).
Linearized Time of Flight (ltof
)¶
Returns a Duration
object corresponding to the linearized time of flight as computed in "Closed loop terminal guidance navigation for a kinetic impactor spacecraft.", Bhaskaran & Kennedy (2014). Acta Astronautica, 103, 322–332 (doi:10.1016/j.actaastro.2014.02.024).^{2}
C3¶
Computes the \(C_3\) of this orbit:
Declination (declination
)¶
Returned in degrees between \([180;180]\).
Distance between two orbits distance_to
¶
The distance in kilometers between two states, computed as follows (where \(r\) and \(r'\) respectively refer to the this
state and the other
state).
Eccentricity (evec
, ecc
)¶
Respectively, the eccentricity vector and its magnitude (i.e. the eccentricity)
For the rest of the MathSpec, let \(e = \mathbf e\).
Eccentric anomaly (ea
)¶
Nyx returns this parameter (\(E\)) in degrees.

Compute the sine of the eccentric anomaly
\[\sin E = \frac{\sqrt{1  e^2} \sin\nu} {1 + e \cos \nu}\] 
Compute the cosine of the eccentric anomaly
\[\cos E = \frac{e + \cos\nu} {1 + e \cos \nu}\] 
Returns the quadrantchecked tangent of those angles (
atan2
), converted to degrees.
Note
This is a conversion from GMAT's StateConversionUtil::TrueToEccentricAnomaly
.
Flight path angle (fpa
)¶

Compute the sine of the FPA, \(f\):
\[ \sin f = \frac{e \sin \nu} {\sqrt{1 + 2e\cos\nu + e^2}} \] 
Compute the cosine of the FPA, \(f\):
\[ \cos f = \frac{1 + e \cos \nu} {\sqrt{1 + 2e\cos\nu + e^2}} \] 
Returns the quadrantchecked tangent of those angles (
atan2
), converted to degrees.
Geodetic height (geodetic_height
)¶
The parameter is returned in kilometers. This is computed using the Vallado approach, Algorithm 12 page 172 in the 4th edition. This accounts for the correction when near the poles. It's a notch complex to write up, so please refer to the code or Vallado for implementation details. As you'll note from the Validation section, it has been validated against Vallado examples.
Warning
This function requires that the orbit already be in a body fixed frame. Nyx will not check that.
Geodetic latitude (geodetic_latitude
)¶
The parameter is returned in degrees between \([180;180]\). This is computed using the Vallado iterative approach, Algorithm 12 page 172 in the 4th edition. It accounts for the flattening of the ellipsoid and its semimajor axis. It's a notch complex to write up, so please refer to the code or Vallado for implementation details. As you'll note from the Validation section, it has been validated against Vallado examples.
Warning
This function requires that the orbit already be in a body fixed frame. Nyx will not check that.
Geodetic longitude (geodetic_longitude
)¶
This is computed using G. Xu and Y. Xu, "GPS", DOI 10.1007/9783662503676_2, 2016, but the validation against the Vallado examples proves to be correct. The following uses the quadrantchecked arctan (atan2
). The angle is returned in degrees and is between \([0;360]\).
Warning
This function requires that the orbit already be in a body fixed frame. Nyx will not check that.
Hyperbolic anomaly (hyperbolic_anomaly
)¶
Computes the hyperbolic anomaly for this hyperbolic orbit in degrees between \([0;360]\), or will return an error (NyxError::NotHyperbolic
) if the orbit is not hyperbolic.
Inclination (inc
)¶
The inclination of this orbit, returned in degrees.
Mean anomaly (ma
)¶

If the eccentricity is below 1, the mean anomaly \(M\) is computed as follows
\[ M = Ee\sin E \] 
If the eccentricity is strictly above 1:
\[ M = \sinh^{1} \frac {\sin\nu \sqrt{e^2  1} }{1 + e \cos\nu} \]
Note
This is a conversion from GMAT's StateConversionUtil::TrueToMeanAnomaly
. This function does support the mean anomaly for hyperbolic orbits (but not for parabolic orbits).
Momentum vector (hvec
, hmag
, hx
, hy
, hz
)¶
The momentum vector, its magnitude, and its X, Y and Z components.
Orbital energy (energy
)¶
The orbit's energy:
Orbital period (period
)¶
The orbital period, returned as highprecision Duration
structure:
Radius of apoapsis (apoapsis
)¶
The radius of apoapsis in kilometers:
Radius of periapsis (periapsis
)¶
The radius of periapsis in kilometers:
Radius unit vector (r_hat
)¶
Unit vector in the direction of the radius vector.
Radius vector (radius
, rmag
)¶
The radius vector, in kilometers.
The magnitude of the radius vector in kilometers.
Right ascension (right_ascension
)¶
Returned in degrees between \([0;360]\).
Right ascension of the ascending node (raan
)¶
The right ascension of the ascending node and is returned in degrees:

Rotate the momentum vector:
\[\begin{equation} \mathbf n = \begin{bmatrix} 0 \\ 0 \\ 1 \\ \end{bmatrix} \times \mathbf h \end{equation}\] 
Compute the RAAN:
\[ \Omega = \cos ^{1} \left( \frac {n_x}{\mathbf n} \right) \] 
Perform a quadrant check: if \(n_y\) is less then zero, return the AoP as \(2\pi  \Omega\) instead.
Semi major axis (sma
)¶
The semi major axis, returned in kilometers.
Semi minor axis (semi_minor_axis
)¶
Returns the semi minor axis in kilometers for hyperbolic and elliptical orbits (will fail for a perfectly circular orbit).
If the eccentricity is less than \(1.0\):
Else:
Semi parameter (semi_parameter
)¶
The semilatus rectum, in kilometers
True anomaly (ta
)¶
The true anomaly is computed as follows and is returned in degrees.

Compute \(\cos \nu\).
We also check that the value is bounded between \([1;1]\) (as it should be mathematically but rounding issues on computers may cause problems): if not, depending on the value of \(\cos \nu\) the phasing is set to either 0 or 180 degrees.
\[ \cos \nu = \mathbf e \cdot \frac {\mathbf r}{e \times \mathbf r } \] 
Compute the true anomaly with a quadrant check.
 If the arccos of \(\nu\) fails (NaN), then a warning is emited and return 0 degrees.
 Else if \(\mathbf r \cdot \mathbf v < 0\), return \(2\pi  \nu\).
 Else, return \(\nu\).
Warning
If the eccentricity is below the eccentricity tolerance, a warning is emitted stating that the true anomaly is illdefined for circular orbits.
True longitude (tlong
)¶
The true longitude is an alternate phasing element. Nyx will return this angle between 0 and 360 degrees.
Velocity declination (velocity_declination
)¶
Returned in degrees between \([180;180]\).
Velocity vector (velocity
, vmag
)¶
The velocity vector, in kilometers per second.
The magnitude of the velocity vector in kilometers per second.
Velocity unit vector (v_hat
)¶
Unit vector in the direction of the velocity vector.
Validation¶
This validation compares the computations of orbital elements in nyx with those in GMAT. Each scenario script is in the subfolder state.
The following table corresponds to the absolute errors between the Nyx computations and those of GMAT. I'll save you the read: the absolute errors are precisely zero for a 64bit floating point representation (double
in C).
Validation
To run all of these test cases, clone the Nyx repo and execute the following command:
cargo test state_def_circ_eq  nocapture
cargo test state_def_circ_inc  nocapture
cargo test state_def_ellip  nocapture
cargo test geodetic  nocapture
Methodology
The validation is trivially the following: a spacecraft is defined in GMAT in the EarthMJ2000Eq frame by its Cartesian state. There is a "Propagate" segment of strictly zero seconds. The desired orbital state computations are exported to a report file.
From a Cartesian state¶
Element / Scenario  circular inclined  circular equatorial  elliptical 

Earth.Energy  0.0  0.0  0.0 
Earth.OrbitPeriod  0.0  0.0  0.0 
Earth.HX  0.0  0.0  0.0 
Earth.HY  0.0  0.0  0.0 
Earth.HZ  0.0  0.0  0.0 
Earth.SMA  0.0  0.0  0.0 
Earth.ECC  0.0  0.0  0.0 
EarthMJ2000Eq.INC  0.0  0.0  0.0 
EarthMJ2000Eq.RAAN  0.0  0.0  0.0 
EarthMJ2000Eq.AOP  0.0  0.0  0.0 
Earth.TA  0.0  0.0  0.0 
Earth.TLONG  0.0  0.0  0.0 
Earth.EA  0.0  0.0  0.0 
Earth.MA  0.0  0.0  0.0 
Earth.RadApo  0.0  0.0  0.0 
Earth.RadPer  0.0  0.0  0.0 
Earth.SemilatusRectum  0.0  0.0  0.0 
From a Keplerian state¶
Element / Scenario  circular inclined  circular equatorial  elliptical 

Earth.X  0.0  0.0  0.0 
Earth.Y  0.0  0.0  0.0 
Earth.Z  0.0  0.0  0.0 
Earth.VX  0.0  0.0  0.0 
Earth.VY  0.0  0.0  0.0 
Earth.VZ  0.0  0.0  0.0 
Earth.Energy  0.0  0.0  0.0 
Earth.OrbitPeriod  0.0  0.0  0.0 
Earth.HX  0.0  0.0  0.0 
Earth.HY  0.0  0.0  0.0 
Earth.HZ  0.0  0.0  0.0 
Earth.SMA  0.0  0.0  0.0 
Earth.ECC  0.0  0.0  0.0 
EarthMJ2000Eq.INC  0.0  0.0  0.0 
EarthMJ2000Eq.RAAN  0.0  0.0  0.0 
EarthMJ2000Eq.AOP  0.0  0.0  0.0 
Earth.TA  0.0  0.0  0.0 
Earth.TLONG  0.0  0.0  0.0 
Earth.EA  0.0  0.0  0.0 
Earth.MA  0.0  0.0  0.0 
Earth.RadApo  0.0  0.0  0.0 
Earth.RadPer  0.0  0.0  0.0 
Earth.SemilatusRectum  0.0  0.0  0.0 
From geodesic¶
Validation
$ cargo test geodetic  nocapture
Compiling nyxspace v1.0.0alpha.1 (/home/chris/Workspace/rust/nyx)
Finished test [unoptimized + debuginfo] target(s) in 35.70s
(...)
Running target/debug/deps/lib48ccd68343ebbebb
running 1 test
[tests/cosmic/state.rs:341] eme2k.semi_major_radius() = 6378.1363
latitude (φ): 0.00e0
longitude (λ): 0.00e0
height: 0.00e0
r_i: 0.00e0
r_j: 0.00e0
r_k: 0.00e0
r_i: 0.00e0
r_j: 0.00e0
r_k: 0.00e0
latitude (φ): 0.00e0
longitude (λ): 0.00e0
height: 0.00e0
test cosmic::state::geodetic_vallado ... ok
Partial (OrbitDual
)¶
Nyx is all about using dual numbers!^{3} An OrbitDual
object can be created from a normal Orbit
structure in order to retrieve the exact partials of many orbital elements with respect to each component of the position and velocity. These can be combined for achieving specific targets and is used for achieving BPlane targets.
List of available partials, always with respect to the position components x,y,z and the velocity components vx, vy, vz. Each of these partials are accessible via their respective wtr_PARAM()
function.
 magnitude of radius vector
rmag
 magnitude of velocity vector
vmag
 X component of the orbital momentum
hx
 Y component of the orbital momentum
hy
 Z component of the orbital momentum
hz
 magnitude of orbital momentum vector
hmag
 orbital energy
energy
 eccentricity
ecc
 inclination
inc
 argument of periapsis
aop
 right ascension of the ascending node
raan
 true anomaly
ta
 true longitude
tlong
 argument of latitude
aol
 periapsis
 apoapsis
 eccentric anomaly
ea
 flight path angle
fpa
 mean anomaly
ma
 semi parameter
semi_parameter
 geodetic longitude
geodetic_longitude
 geodetic latitude
geodetic_latitude
 geodetic height
geodetic_height
 right ascension
right_ascension
 decliation
decliation
 semi minor axis
semi_minor_axis
 velocity decliation
velocity_decliation
 \(C_3\)
c3
 hyperbolic anomaly
hyperbolic_anomaly

Nyx allows initialization from geodesic elements only for the following celestial bodies: Mercury, Venus, Earth, Luna/Earth Moon, and Mars. Although the Jupiter, Saturn, Uranus, and Neptune also have an angular velocity defined in Nyx, they do not have an ellipsoid flatenning parameter. ↩

Two other computation methods were attempted (McMahon and Jah). The first prevented convergence of the shooting algorithm and the second led to near infinite LTOF. ↩

Please refer to Dual Numbers for a primer on dual number theory. ↩