orbitalMotion

Defines

N_DEBYE_PARAMETERS

jPerturb(...)

Enums

enum CelestialObject_t

Values:

enumerator CELESTIAL_MERCURY

enumerator CELESTIAL_VENUS

enumerator CELESTIAL_EARTH

enumerator CELESTIAL_MOON

enumerator CELESTIAL_MARS

enumerator CELESTIAL_PHOBOS

enumerator CELESTIAL_DEIMOS

enumerator CELESTIAL_JUPITER

enumerator CELESTIAL_SATURN

enumerator CELESTIAL_URANUS

enumerator CELESTIAL_NEPTUNE

enumerator CELESTIAL_PLUTO

enumerator CELESTIAL_SUN

enumerator MAX_CELESTIAL

Functions

double E2f(double E, double e): Purpose: Maps eccentric anomaly angles into true anomaly angles. This function requires the orbit to be either circular or non-rectilinar elliptic orbit. Inputs: Ecc = eccentric anomaly (rad) e = eccentricity (0 <= e < 1) Outputs: f = true anomaly (rad)

double E2M(double E, double e): Purpose: Maps the eccentric anomaly angle into the corresponding mean elliptic anomaly angle. Both 2D and 1D elliptic orbit are allowed. Inputs: Ecc = eccentric anomaly (rad) e = eccentricity (0 <= e < 1) Outputs: M = mean elliptic anomaly (rad)

double f2E(double f, double e): Purpose: Maps true anomaly angles into eccentric anomaly angles. This function requires the orbit to be either circular or non-rectilinar elliptic orbit. Inputs: f = true anomaly angle (rad) e = eccentricity (0 <= e < 1) Outputs: Ecc = eccentric anomaly (rad)

double f2H(double f, double e): Purpose: Maps true anomaly angles into hyperbolic anomaly angles. This function requires the orbit to be hyperbolic Inputs: f = true anomaly angle (rad) e = eccentricity (e > 1) Outputs: H = hyperbolic anomaly (rad)

double H2f(double H, double e): Purpose: Maps hyperbolic anomaly angles into true anomaly angles. This function requires the orbit to be hyperbolic Inputs: H = hyperbolic anomaly (rad) e = eccentricity (e > 1) Outputs: f = true anomaly angle (rad)

double H2N(double H, double e): Purpose: Maps the hyperbolic anomaly angle H into the corresponding mean hyperbolic anomaly angle N. Inputs: H = hyperbolic anomaly (rad) e = eccentricity (e > 1) Outputs: N = mean hyperbolic anomaly (rad)

double M2E(double M, double e)

Solve the elliptic Kepler equation using bracketed Newton iteration.

Note

The result preserves the revolutions in M. Very large unwrapped anomalies are limited by input/output angle resolution.

Parameters:

M – Finite mean elliptic anomaly [rad].
e – Finite eccentricity, 0 <= e < 1.

Returns:

Eccentric anomaly [rad], or NaN with an error diagnostic for invalid inputs or failure to converge within 200 iterations.

double N2H(double N, double e)

Solve the hyperbolic Kepler equation using bracketed Newton iteration.

Note

The equation is scaled to avoid overflow for large finite mean anomalies and evaluated without subtractive cancellation near the parabolic limit.

Parameters:

N – Finite mean hyperbolic anomaly [rad].
e – Finite eccentricity, e > 1.

Returns:

Hyperbolic anomaly [rad], or NaN with an error diagnostic for invalid inputs or failure to converge within 200 iterations.

void elem2rv(double mu, ClassicElements *elements, double *rVec, double *vVec)

Convert classical elements to a finite inertial Cartesian state.

Note

Inclination must satisfy 0 <= i <= pi. Ellipses use a > 0, e < 1; hyperbolas use a < 0, e > 1. No eccentricity tolerance selects a conic. Planar parabolas use e = 1, a = 0, rPeriap > 0 and true anomaly f.

Note

Rectilinear ellipses use e = 1, a > 0 and f = eccentric anomaly E, with outward motion for sin(E) > 0. Rectilinear hyperbolas use e = 1, a < 0 and f = signed hyperbolic anomaly H, outward for H > 0. Rectilinear parabolas use e = 1, a = rPeriap = 0, rmag > 0 and f = 0 for outward motion or f = pi for inward motion. For radial states the perifocal x axis is the outward radial ray; the plane about it is arbitrary. Extra fields are read only where required by these conventions.

Note

Invalid/nonfinite inputs, collisions, open-conic anomalies at or beyond infinity, and numerical overflow produce a diagnostic and all-NaN vectors.

Parameters:

mu – Positive gravitational parameter [L^3/s^2].
elements – Orbital elements, with distances in L and angles in radians.
rVec – [out] Inertial position vector [L].
vVec – [out] Inertial velocity vector [L/s].

void rv2elem(double mu, double *rVec, double *vVec, ClassicElements *elements)

Convert a Cartesian state to classical orbital elements.

Note

Inputs must be finite and use consistent distance units. Inclination is in [0, pi]; Omega, omega, and nonradial f are in [0, 2*pi), including hyperbolic f. Degenerate conics use the conventions documented in elem2rv().

Note

A dimensionless tolerance of 1e-12 selects singular-angle conventions. For e <= tolerance, omega = 0 and f is the argument of latitude. For sin(i) <= tolerance, Omega = 0 and i is set to 0 or pi; omega measures periapsis from the inertial x axis in the direction of orbital motion. If both conditions hold, f measures position from that axis in the direction of motion. Replacing unresolved directions introduces relative state errors on the order of the singularity tolerance.

Note

Parabolic energy and zero transverse velocity are detected within 8 machine epsilons of the local energy and speed scales, respectively. Radial states use a deterministic plane about the outward radial ray. Near zero energy, relative accuracy of a and state reconstruction from rounded a, e degrades as machine epsilon divided by abs(1-e). Nonradial finite-energy states whose eccentricity rounds to one or the wrong side of one cannot be encoded consistently. These, invalid inputs and numerical overflow produce a diagnostic and all-NaN elements.

Parameters:

mu – Positive gravitational parameter [L^3/s^2].
rVec – Nonzero inertial position vector [L].
vVec – Inertial velocity vector [L/s].
elements – [out] Orbital elements, with distances in L and angles in radians.

void clMeanOscMap(double req, double J2, ClassicElements *elements, ClassicElements *elements_p, double sgn)

Apply the first-order J2 mean/osculating element mapping.

Note

See clMeanOscMapChecked() for the supported domain. Rejected domains retain historical results. Checked calls provide explicit domain validation. Invalid inclinations are always rejected; use 0 <= i <= pi.

Parameters:

req – Equatorial radius, in the same length unit as elements->a.
J2 – Dimensionless second zonal harmonic coefficient.
elements – Input classical elements, with angles in radians.
elements_p – Output classical elements, with the same units as the input.
sgn – +1 for mean to osculating, -1 for osculating to mean.

void clMeanOscMapChecked(double req, double J2, ClassicElements *elements, ClassicElements *elements_p, double sgn)

Apply the first-order J2 mean/osculating element mapping.

Note

Follows Appendix F of Analytical Mechanics of Space Systems, Schaub and Junkins, fourth edition. Inputs must be finite and elliptic with inclination in [0, pi], away from the critical inclinations where 1 - 5*cos(i)^2 = 0. The perturbation must be small; opposite mappings leave second-order J2 errors.

Note

Both equatorial limits are supported. Whole revolutions in true anomaly do not affect the mapped Cartesian state; output angles need not retain input revolution counts. J2 = 0 copies the six classical elements exactly. Only a, e, i, Omega, omega and f are written; auxiliary fields are unchanged. Input and output may alias.

Note

Requires req > 0 and sgn = +/-1. For nonzero J2, rejects |1 - 5*cos(i)^2| <= sqrt(DBL_EPSILON) and nonphysical or nonfinite mapped elements. This numerical guard does not establish perturbation-model accuracy. Invalid inputs or results produce a diagnostic and six NaN classical elements.

Parameters:

req – Equatorial radius, in the same length unit as elements->a.
J2 – Dimensionless second zonal harmonic coefficient.
elements – Input classical elements, with angles in radians.
elements_p – Output classical elements, with the same units as the input.
sgn – +1 for mean to osculating, -1 for osculating to mean.

void clElem2eqElem(ClassicElements *elements_cl, equinoctialElements *elements_eq)

Map finite elliptic classical elements to prograde equinoctial elements.

Note

See clElem2eqElemChecked() for the supported domain. Rejected domains retain historical results. Checked calls provide explicit domain validation. Invalid inclinations are always rejected; use 0 <= i <= pi.

Parameters:

elements_cl – Input with a > 0, 0 <= e < 1 and 0 <= i < pi; angles in radians.
elements_eq – Output with the same length unit as the input and angles in radians.

void clElem2eqElemChecked(ClassicElements *elements_cl, equinoctialElements *elements_eq)

Map finite elliptic classical elements to prograde equinoctial elements.

Note

The retrograde equator is singular. Rejects cos(i/2) <= 8*DBL_EPSILON. Invalid inputs or nonfinite results produce a diagnostic and seven NaN outputs.

Parameters:

elements_cl – Input with a > 0, 0 <= e < 1 and 0 <= i < pi; angles in radians.
elements_eq – Output with the same length unit as the input and angles in radians.

void hillFrame(double *rc_N, double *vc_N, double HN[3][3])

Compute the Hill DCM from a nonrectilinear chief state.

Note

See hillFrameChecked() for the supported domain. Rejected domains retain historical results. Checked calls provide explicit domain validation.

Parameters:

rc_N – Finite inertial position [L], with finite nonzero norm.
vc_N – Finite inertial velocity [L/s], with finite nonzero norm.
HN – Output DCM mapping inertial components to Hill components.

void hillFrameChecked(double *rc_N, double *vc_N, double HN[3][3])

Compute the Hill DCM from a nonrectilinear chief state.

Note

Requires norm(cross(unit(rc_N), unit(vc_N))) > 8*DBL_EPSILON. Invalid states produce a diagnostic and an all-NaN matrix.

Parameters:

rc_N – Finite inertial position [L], with finite nonzero norm.
vc_N – Finite inertial velocity [L/s], with finite nonzero norm.
HN – Output DCM mapping inertial components to Hill components.

void hill2rv(double *rc_N, double *vc_N, double *rho_H, double *rhoPrime_H, double *rd_N, double *vd_N)

Map a deputy state from the rotating chief Hill frame to inertial coordinates.

Note

See hill2rvChecked() for the supported domain. Rejected domains retain historical results. Checked calls provide explicit domain validation.

Parameters:

rc_N – Chief inertial position [L], satisfying hillFrameChecked() requirements.
vc_N – Chief inertial velocity [L/s], satisfying hillFrameChecked() requirements.
rho_H – Finite Hill relative position [L].
rhoPrime_H – Finite rotating-frame relative velocity [L/s].
rd_N – Output deputy inertial position [L].
vd_N – Output deputy inertial velocity [L/s].

void hill2rvChecked(double *rc_N, double *vc_N, double *rho_H, double *rhoPrime_H, double *rd_N, double *vd_N)

Map a deputy state from the rotating chief Hill frame to inertial coordinates.

Note

Invalid states, unrepresentable positive angular rate, or nonfinite outputs produce a diagnostic and two all-NaN vectors. Use consistent distance units.

Parameters:

rc_N – Chief inertial position [L], satisfying hillFrame() requirements.
vc_N – Chief inertial velocity [L/s], satisfying hillFrame() requirements.
rho_H – Finite Hill relative position [L].
rhoPrime_H – Finite rotating-frame relative velocity [L/s].
rd_N – Output deputy inertial position [L].
vd_N – Output deputy inertial velocity [L/s].

void rv2hill(double *rc_N, double *vc_N, double *rd_N, double *vd_N, double *rho_H, double *rhoPrime_H)

Express the deputy state relative to the chief in its rotating Hill frame.

Note

See rv2hillChecked() for the supported domain. Rejected domains retain historical results. Checked calls provide explicit domain validation.

Parameters:

rc_N – Chief inertial position [L], satisfying hillFrameChecked() requirements.
vc_N – Chief inertial velocity [L/s], satisfying hillFrameChecked() requirements.
rd_N – Finite deputy inertial position [L].
vd_N – Finite deputy inertial velocity [L/s].
rho_H – Output Hill relative position [L].
rhoPrime_H – Output rotating-frame relative velocity [L/s].

void rv2hillChecked(double *rc_N, double *vc_N, double *rd_N, double *vd_N, double *rho_H, double *rhoPrime_H)

Express the deputy state relative to the chief in its rotating Hill frame.

Note

Invalid states, unrepresentable positive angular rate, or nonfinite outputs produce a diagnostic and two all-NaN vectors. Use consistent distance units.

Parameters:

rc_N – Chief inertial position [L], satisfying hillFrame() requirements.
vc_N – Chief inertial velocity [L/s], satisfying hillFrame() requirements.
rd_N – Finite deputy inertial position [L].
vd_N – Finite deputy inertial velocity [L/s].
rho_H – Output Hill relative position [L].
rhoPrime_H – Output rotating-frame relative velocity [L/s].

double atmosphericDensity(double alt)

Purpose: This program computes the atmospheric density based on altitude supplied by user. This function uses a curve fit based on atmospheric data from the Standard Atmosphere 1976 Data. This function is valid for altitudes ranging from 100km to 1000km.

Note: This code can only be applied to spacecraft orbiting the Earth Inputs: alt = altitude in km Outputs: density = density at the given altitude in kg/m^3

double debyeLength(double alt)

Compute the Debye length for a given altitude.

Note

Values above 1000 km are highly speculative.

Parameters:: alt – Altitude [km], valid from 200 km through 35000 km.
Returns:: Debye length [m], or NaN if alt is non-finite or outside the valid range.

void atmosphericDrag(double Cd, double A, double m, double *rvec, double *vvec, double *advec)

Estimate Earth’s atmospheric drag acceleration opposite the supplied velocity.

Note

Assumes an atmosphere at rest in the input inertial frame; no rotation or wind correction is applied. Uses atmosphericDensity(), whose fit is valid from 100 to 1000 km; positive altitudes below this range are extrapolated, and above 1000 km the existing exponential tail is retained, without a hard cutoff.

Note

Zero speed, area, coefficient, or computed density gives zero acceleration after input validation. Invalid inputs, unrepresentable norms, or acceleration overflow produce a diagnostic and three NaNs. Acceleration underflow is allowed. The output may alias either input vector; all arrays must hold three doubles.

Parameters:

Cd – Finite, nonnegative drag coefficient [-].
A – Finite, nonnegative projected area [m^2].
m – Finite, positive spacecraft mass [kg].
rvec – Finite inertial position three-vector [km], above Earth’s equatorial radius.
vvec – Finite inertial velocity three-vector [km/s].
advec – [out] Drag acceleration in the input frame [km/s^2].

void jPerturb(double *rvec, int num, double *ajtot, ...)

Retain the legacy variadic symbol for explicit-planet native callers.

Note

Exactly one trailing int planet identifier is required. Header-based three-argument calls select Earth through the header adapter without reading varargs.

Parameters:

rvec – Position [km].
num – Maximum zonal degree.
ajtot – Output acceleration [km/s^2].

void jPerturbEarth(double *rvec, int num, double *ajtot)

Compute Earth’s zonal perturbation acceleration through the requested degree.

Note

This fixed three-argument entry point selects Earth. For another body, use jPerturbForPlanet(rvec, num, ajtot, planetID). See jPerturbForPlanet() for model assumptions and failure behavior.

Parameters:

rvec – Finite, nonzero planet-centered position [km], with z along the spin axis.
num – Maximum included degree, from 2 through 6 (inclusive).
ajtot – [out] Perturbation acceleration in the input frame [km/s^2].

void jPerturbLegacy(double *rvec, int num, double *ajtot, int planetID)

Explicit-planet adapter preserving the original jPerturb conventions.

Note

Sun/Pluto retain the historical unchanged-output behavior. Use jPerturbForPlanet() for checked selection and explicit zero-zonal outputs.

Parameters:

rvec – Position [km].
num – Maximum zonal degree.
ajtot – Output acceleration [km/s^2].
planetID – Legacy planet identifier; unknown values select Earth.

void jPerturbForPlanet(double *rvec, int num, double *ajtot, CelestialObject_t planetID)

Compute a selected body’s zonal perturbation acceleration.

Note

Earth includes J2 through J6. Mercury, Venus, Moon, Mars, Jupiter, Saturn, Uranus and Neptune include J2 only. Sun and Pluto use a zero-zonal approximation and return a zero vector. Phobos and Deimos are unsupported. These approximations omit central gravity, tesseral terms and other forces.

Note

The external gravity expansion is intended for positions outside the body’s reference radius; the surface boundary is not enforced. All inputs are validated even for zero-zonal models. Invalid inputs or nonfinite results produce a diagnostic and an all-NaN vector. Input/output arrays may alias.

Parameters:

rvec – Finite, nonzero planet-centered position [km], with z along the spin axis.
num – Maximum included degree, from 2 through 6 (inclusive).
ajtot – [out] Perturbation acceleration in the input frame [km/s^2].
planetID – Body whose coefficients are used; unsupported identifiers are errors.

void solarRad(double A, double m, double *sunvec, double *arvec)

Purpose: Computes the inertial solar radiation force vectors based on cross-sectional Area and mass of the spacecraft and the position vector of the planet to the sun. Note: It is assumed that the solar radiation pressure decreases quadratically with distance from sun (in AU)

Solar Radiation Equations obtained from Earth Space and Planets Journal Vol. 51, 1999 pp. 979-986 Inputs: A = Cross-sectional area of the spacecraft that is facing the sun in m^2. m = The mass of the spacecraft in kg. sunvec = Position vector to the Sun in units of AU. Earth has a distance of 1 AU. Outputs: arvec = The inertial acceleration vector due to the effects of Solar Radiation pressure in km/sec^2. The vector components of the output are the same as the vector components of the sunvec input vector.

static inline void bskJPerturbDispatch(double *rvec, int num, double *ajtot, int planetID, ...)

Header adapter for three- and four-argument jPerturb source calls.

Note

Trailing arguments are ignored. Forwarding the complete argument list avoids counting commas inside C compound literals or C++ braced initializers.

Parameters:

rvec – Position [km].
num – Maximum zonal degree.
ajtot – Output acceleration [km/s^2].
planetID – Explicit planet or the Earth default appended by the macro.

struct ClassicElements

#include <orbitalMotion.h>

Classical orbit elements, using consistent distance units L and angles in radians.

Note

See elem2rv() for planar parabolic and rectilinear encodings. In particular, e = 1 with nonzero a denotes radial motion, not a planar parabola.

Public Members

double a: object semi-major axis

double e: Eccentricity of the orbit.

double i: inclination of the orbital plane

double Omega: Right ascension of the ascending node.

double omega: Argument of periapsis of the orbit.

double f: True anomaly; radial E/H or parabolic direction flag as documented in elem2rv().

double rmag: Position magnitude [L]; required input for a rectilinear parabola.

double alpha: Inverted semi-major axis (extra).

double rPeriap: Periapsis radius [L]; required input when a = 0, e = 1.

double rApoap: Apoapsis radius [L]; signed algebraic value for hyperbolas, zero for parabolas.

struct equinoctialElements

#include <orbitalMotion.h>

This structure contains the set of Keplerian orbital elements that define the spacecraft translational state. It is operated on by the orbital element routines and the OrbElemConvert module.

equinoctial elment struct definition

Public Members

double a: semi-major axis

double P1: e*sin(omega+Omega)

double P2: e*cos(omega+Omega)

double Q1: tan(i/2)*sin(Omega)

double Q2: tan(i/2)*cos(Omega)

double l: Omega+omega+M.

double L: Omega+omega+f.