orbitalMotion

orbitalMotion.E2M(Ecc: float, e: float) → float

Maps the eccentric anomaly angle into the corresponding mean elliptic anomaly angle. Both 2D and 1D elliptic orbit are allowed.

Parameters:

• Ecc – eccentric anomaly (rad)

• e – eccentricity (0 <= e < 1)

Returns:

M, mean elliptic anomaly (rad)

orbitalMotion.E2f(Ecc: float, e: float) → float

Maps eccentric anomaly angles into true anomaly angles This function requires the orbit to be either circular or non-rectilinear elliptic orbit

Parameters:

• Ecc – eccentric anomaly (rad)

• e – eccentric (0 <= e < 1)

Returns:

f, true anomaly (rad)

orbitalMotion.H2N(H: float, e: float) → float

Maps the hyperbolic anomaly angle H into the corresponding mean hyperbolic anomaly angle N.

Parameters:

• H – hyperbolic anomaly (rad)

• e – eccentricity (e > 1)

Returns:

N, mean hyperbolic anomaly (rad)

orbitalMotion.H2f(H: float, e: float) → float

Maps hyperbolic anomaly angles into true anomaly angles. This function requires the orbit to be hyperbolic

Parameters:

• H – hyperbolic anomaly (rad)

• e – eccentricity (e > 1)

Returns:

f, true anomaly angle (rad)

orbitalMotion.M2E(M: float, e: float) → float

Solve the elliptic Kepler equation using bracketed Newton iteration.

The returned eccentric anomaly preserves the revolutions in M. Very large unwrapped anomalies are limited by the resolution of the input and output angles.

Parameters:

• M – mean elliptic anomaly (rad)

• e – eccentricity (0 <= e < 1)

Returns:

eccentric anomaly (rad); NumPy scalar types and zero-dimensional mean-array results are preserved

Raises:

• ValueError – if either input is non-finite or eccentricity is outside its domain

• RuntimeError – if the iteration limit is reached without convergence

orbitalMotion.N2H(N: float, e: float) → float

Solve the hyperbolic Kepler equation using bracketed Newton iteration.

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

Parameters:

• N – mean hyperbolic anomaly (rad)

• e – eccentricity (e > 1)

Returns:

hyperbolic anomaly (rad); NumPy scalar types and zero-dimensional mean-array results are preserved

Raises:

• ValueError – if either input is non-finite or eccentricity is outside its domain

• RuntimeError – if the iteration limit is reached without convergence

orbitalMotion.atmosphericDensity(alt: float) → float

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

Parameters:

alt – altitude in km

Returns:

density at the given altitude in kg/m^3

orbitalMotion.atmosphericDrag(Cd: float, A: float, m: float, rvec: ndarray, vvec: ndarray) → ndarray

Estimate Earth’s drag acceleration opposite the supplied velocity.

The atmosphere is assumed at rest in the input inertial frame; no rotation or wind correction is applied. The atmosphericDensity() fit is valid from 100 to 1000 km. Positive altitudes below this range are extrapolated; above 1000 km its existing exponential tail is retained without a cutoff.

Zero speed, area, coefficient, or computed density gives a zero vector after input validation. Acceleration underflow is allowed. Inputs are not modified; sequences and NumPy arrays with shape (3,) are accepted.

Parameters:

• Cd – finite, nonnegative drag coefficient (dimensionless)

• A – finite, nonnegative projected area (m^2)

• m – finite, positive spacecraft mass (kg)

• rvec – finite inertial position three-vector above Earth’s equatorial radius (km)

• vvec – finite inertial velocity three-vector (km/s)

Returns:

drag acceleration in the input frame (km/s^2), as a NumPy three-vector

Raises:

ValueError – for invalid inputs, unrepresentable norms, or acceleration overflow

orbitalMotion.clElem2eqElem(elements_cl: ClassicElements, elements_eq: EquinoctialElements)

Convert elements with the compatibility policy described in clMeanOscMap().

Parameters and supported-domain results are those of clElem2eqElemChecked().

orbitalMotion.clElem2eqElemChecked(elements_cl: ClassicElements, elements_eq: EquinoctialElements)

Convert an ellipse to prograde equinoctial elements.

Parameters:

• elements_cl – finite classical elements with a > 0, 0 <= e < 1, 0 <= i < pi; angles in radians

• elements_eq – output elements, using the input length unit and radians

Returns:

None; writes a, P1, P2, Q1, Q2, l, L to elements_eq

Raises:

ValueError – for invalid inputs, the retrograde-equatorial singularity, or nonfinite output; the output object is unchanged on failure

The prograde equator is supported. The retrograde equator is singular: cos(i/2) <= 8*machine epsilon is rejected. Close to this boundary, the coordinates remain ill-conditioned.

orbitalMotion.clMeanOscMap(req: float, J2: float, oe: ClassicElements, oep: ClassicElements, sign: int)

Map mean/osculating elements while preserving historical failure behavior.

Parameters and supported-domain results are those of clMeanOscMapChecked(). Inputs rejected by that function use the original formulas, including their historical exceptions or nonfinite results. Use the checked function to opt into explicit domain validation. Inclination must be finite and in [0, pi], including on the legacy path.

orbitalMotion.clMeanOscMapChecked(req: float, J2: float, oe: ClassicElements, oep: ClassicElements, sign: int)

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

The mapping 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], and sufficiently far from the critical inclinations where 1 - 5*cos(i)**2 = 0. The perturbation must be small; applying the opposite mapping leaves errors of second order in J2.

Both equatorial limits are supported. Whole revolutions in true anomaly do not change the mapped Cartesian state; output angles need not retain input revolution counts. With J2 = 0, the six classical elements are copied exactly, including their angle branches. Only a, e, i, Omega, omega, f are written; auxiliary fields are unchanged. Input and output may be the same object.

Parameters:

• req – equatorial radius, in the same length unit as oe.a

• J2 – dimensionless second zonal harmonic coefficient

• oe – input classical elements, with angles in radians

• oep – output classical elements, with the same units as oe

• sign – +1 for mean to osculating, -1 for osculating to mean

Returns:

None; the mapped elements are written to oep

Raises:

ValueError – for invalid inputs, singular inclinations or invalid outputs

For nonzero J2, abs(1 - 5*cos(i)**2) <= sqrt(machine epsilon) is rejected. Mapped elements must remain finite and elliptic. This numerical guard does not establish the accuracy of the perturbation model. On failure, oep is unchanged. Zero J2 bypasses the critical-inclination guard only.

orbitalMotion.debyeLength(alt: float) → float

This program computes the debyeLength length for a given altitude and is valid for altitudes ranging from 200 km to GEO (35000km). However, all values above 1000 km are HIGHLY speculative at this point.

Parameters:

alt – altitude in km

Returns:

debye length given in m

orbitalMotion.elem2rv(mu: float, elements: ClassicElements) → Tuple[ndarray, ndarray]

Convert classical elements to a finite inertial Cartesian state.

Inclination must be in [0, pi]. Distances use the same unit L as mu. Ellipses have a > 0, e < 1; hyperbolas have a < 0, e > 1. No eccentricity tolerance selects a conic. The C and Python core conversions share these degenerate conventions:

• Planar parabola: e = 1, a = 0, rPeriap > 0; f is true anomaly.

• Rectilinear ellipse: e = 1, a > 0; f is eccentric anomaly E, with outward motion for sin(E) > 0.

• Rectilinear hyperbola: e = 1, a < 0; f is signed hyperbolic anomaly H, with outward motion for H > 0.

• Rectilinear parabola: e = 1, a = rPeriap = 0; rmag > 0 supplies radius, and f = 0 selects outward motion or f = pi inward motion.

For rectilinear states, i, Omega, omega orient the outward radial ray along the perifocal x axis; the plane about that ray is arbitrary. Collisions and open-conic anomalies at or beyond infinity are rejected. The extra fields are read only where required by these conventions.

Parameters:

• mu – positive gravitational parameter (L^3/s^2)

• elements – orbital elements, with distances in L and angles in radians

Returns:

inertial position (L) and velocity (L/s)

Raises:

ValueError – if the elements do not describe a finite physical state

orbitalMotion.elem2rv_parab(mu: float, elements: ClassicElements) → Tuple[ndarray, ndarray]

Convert elements using the legacy Python parabolic encoding.

This adapter accepts e = 1, a < 0 for a planar parabola with periapsis radius -a, independently of auxiliary fields such as rPeriap. Use elem2rv() for rectilinear hyperbolas. Other conventions follow it, except for rectilinear elliptic motion. This adapter retains inward velocity for positive sin(E) on rectilinear ellipses. This convention remains supported. The input object is not modified.

Parameters:

• mu – positive gravitational parameter (L^3/s^2)

• elements – orbital elements, with distances in L and angles in radians

Returns:

inertial position (L) and velocity (L/s)

Raises:

ValueError – if the elements do not describe a finite physical state

orbitalMotion.f2E(f: float, e: float) → float

Maps true anomaly angles into eccentric anomaly angles. This function requires the orbit to be either circular or non-rectilinear elliptic orbit.

Parameters:

• f – true anomaly angle (rad)

• e – eccentricity (0 <= e < 1)

Returns:

Ecc, eccentric anomaly (rad)

orbitalMotion.f2H(f: float, e: float) → float

Maps true anomaly angles into hyperbolic anomaly angles. This function requires the orbit to be hyperbolic

Parameters:

• f – true anomaly angle (rad)

• e – eccentricity (e > 1)

Returns:

H, hyperbolic anomaly (rad)

orbitalMotion.hill2rv(rc_N: ndarray, vc_N: ndarray, rho_H: ndarray, rhoPrime_H: ndarray) → Tuple[ndarray, ndarray]

Convert from Hill coordinates with the compatibility policy of hillFrame().

Parameters and supported-domain results are those of hill2rvChecked().

orbitalMotion.hill2rvChecked(rc_N: ndarray, vc_N: ndarray, rho_H: ndarray, rhoPrime_H: ndarray) → Tuple[ndarray, ndarray]

Map a deputy state in the chief’s rotating Hill frame to inertial vectors.

Parameters:

• rc_N – chief inertial position three-vector, in a consistent length unit

• vc_N – chief inertial velocity three-vector, in length units per second

• rho_H – Hill relative position three-vector, in the same length unit

• rhoPrime_H – Hill rotating-frame velocity three-vector, in the same velocity unit

Returns:

deputy inertial position and velocity vectors

Raises:

ValueError – for invalid vectors, an unresolved chief plane (see hillFrame()), or unrepresentable angular rate or outputs

orbitalMotion.hillFrame(rc_N: ndarray, vc_N: ndarray) → ndarray

Compute a Hill frame while preserving historical behavior for invalid states.

Parameters and supported-domain results are those of hillFrameChecked(). Use the checked function to opt into explicit domain validation.

orbitalMotion.hillFrameChecked(rc_N: ndarray, vc_N: ndarray) → ndarray

Compute the Hill frame DCM HN for a nonrectilinear chief state.

Parameters:

• rc_N – finite inertial position three-vector, in a consistent length unit

• vc_N – finite inertial velocity three-vector, in length units per second

Returns:

DCM mapping inertial frame N components to Hill frame H components

Raises:

ValueError – for invalid or zero vectors, unrepresentable norms, or norm(cross(unit(rc_N), unit(vc_N))) <= 8*machine epsilon

orbitalMotion.jPerturb(rvec: ndarray, num: int, planet: str = 'CELESTIAL_EARTH', *, strict: bool = False) → ndarray | None

Compute a body’s zonal perturbation acceleration through degree num.

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 with strict=True. Phobos, Deimos and unknown body names are unsupported by strict=True. The default retains the original Earth fallback, None for Sun/Pluto, and fractional-degree truncation.

This external gravity expansion is intended for positions outside the body’s reference radius; the surface boundary is not enforced. The input frame must be planet-centered, with its z axis along the spin axis. Central gravity, tesseral terms and other forces are not included. With strict=True, all inputs are validated, including zero-zonal models.

Parameters:

• rvec – finite, nonzero position vector with shape (3,), in kilometers

• num – integer maximum included degree, from 2 through 6 (inclusive)

• planet – body name prefixed by CELESTIAL_; defaults to Earth

• strict – opt in to checked body/degree selection and zero-vector Sun/Pluto results; the default preserves the original dispatch behavior

Returns:

perturbation acceleration in the input frame, in km/s^2, or the legacy None sentinel for Sun/Pluto when strict=False

Raises:

ValueError – for invalid inputs, unsupported bodies or nonfinite results

orbitalMotion.rv2elem(mu: float, rVec: ndarray, vVec: ndarray, *, canonical: bool = False) → ClassicElements

Convert a Cartesian state to classical orbital elements.

Inputs must be finite and use consistent distance units. Inclination i is in [0, pi]; Omega and omega are in [0, 2*pi). Nonradial true anomaly f is in [0, 2*pi) for ellipses/parabolas and [-pi, pi] for hyperbolas. With canonical=True, degenerate conics use the conventions of elem2rv(). The default retains the supported Python parabolic encoding a = -rPeriap, rApoap = -1.

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. Equatorial retrograde longitudes therefore have the opposite sign to inertial azimuths.

• 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 this tolerance. Energy is classified as parabolic only within 8 machine epsilons of the local kinetic/potential energy scale. Radial motion is detected within 8 machine epsilons of zero transverse velocity; the plane is chosen deterministically about the outward 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 and are rejected.

Parameters:

• mu – positive gravitational parameter (L^3/s^2)

• rVec – nonzero inertial position vector (L)

• vVec – inertial velocity vector (L/s)

• canonical – if True, use the C-compatible parabolic encoding; otherwise retain the original Python a = -rPeriap, rApoap = -1 result

Returns:

orbital elements, with distances in L and angles in radians; a, e, rmag, alpha, rPeriap and rApoap retain NumPy float64 scalar arithmetic

Raises:

ValueError – if the state is invalid or its elements are unrepresentable

orbitalMotion.rv2elem_parab(mu: float, rVec: ndarray, vVec: ndarray) → ClassicElements

Convert a Cartesian state using the legacy Python parabolic encoding.

The computation and singular-angle conventions are those of rv2elem(). Planar parabolas return a = rApoap = -rPeriap. Rectilinear states have rPeriap = 0; use elem2rv() to reconstruct radial hyperbolas because elem2rv_parab() retains its negative-a planar encoding. Node, periapsis and nonradial true anomaly retain signed angles in [-pi, pi]; circular orbits place periapsis at the current position. Radial elliptic anomaly retains the supported inward-positive convention.

Parameters:

• mu – positive gravitational parameter (L^3/s^2)

• rVec – nonzero inertial position vector (L)

• vVec – inertial velocity vector (L/s)

Returns:

orbital elements, with distances in L and angles in radians

Raises:

ValueError – if the state is invalid or its elements are unrepresentable

orbitalMotion.rv2hill(rc_N: ndarray, vc_N: ndarray, rd_N: ndarray, vd_N: ndarray) → Tuple[ndarray, ndarray]

Convert to Hill coordinates with the compatibility policy of hillFrame().

Parameters and supported-domain results are those of rv2hillChecked().

orbitalMotion.rv2hillChecked(rc_N: ndarray, vc_N: ndarray, rd_N: ndarray, vd_N: ndarray) → Tuple[ndarray, ndarray]

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

Parameters:

• rc_N – chief inertial position three-vector, in a consistent length unit

• vc_N – chief inertial velocity three-vector, in length units per second

• rd_N – deputy inertial position three-vector, in the same length unit

• vd_N – deputy inertial velocity three-vector, in the same velocity unit

Returns:

Hill relative position and rotating-frame velocity vectors

Raises:

ValueError – for invalid vectors, an unresolved chief plane (see hillFrame()), or unrepresentable angular rate or outputs

orbitalMotion.solarRad(A: float, m: float, sunvec: ndarray) → ndarray

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

Parameters:

• 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.

Returns:

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.