Orbit Propagation Methods

The Aerospace Toolbox supports two top-level orbit propagation methods: Kepler (unperturbed) and Numerical (high precision).

Kepler (unperturbed)

Kepler orbit propagation is the process of numerically computing and predicting the position and velocity of an object in space, based on Kepler's laws of planetary motion. Kepler's laws provide a mathematical description of the motion of objects in elliptical orbits around a central body, such as a planet orbiting the Sun.

To propagate a Keplerian orbit, you typically need these elements:

orbital elements

In this diagram, the orbital plane (yellow) intersects a reference plane (gray). For Earth-orbiting satellites, the reference plane is usually the IJ-plane of the GCRF.

These two elements define the shape and size of the ellipse:

Eccentricity (e) — Describe how elongated the shape of the ellipse is when compared to a circle.
Semimajor axis (a) — The sum of the periapsis and apoapsis distances divided by 2. Periapsis is the point at which an orbiting object is closest to the center of mass of the body it is orbiting. Apoapsis is the point at which an orbiting object is farthest away from the center of mass of the body it is orbiting. For classic two-body orbits, the semimajor axis is the distance between the centers of the bodies.

These two elements define the orientation of the orbital plane in which the ellipse is embedded:

Inclination (i) — The vertical tilt of the ellipse with respect to the reference plane measured at the ascending node. The ascending node is the location where the orbit passes upward through the reference plane (the green angle i in the diagram). The tilt angle is measured perpendicular to the line of intersection between orbital plane and the reference plane. Any three points on an ellipse defines the ellipse orbital plane.
Starting with an equatorial orbit, the orbital plane can be tilted up. The angle it is tilted up from the equator is referred to as the inclination angle, i, in the range [0,180]. Because the center of the Earth must always be in the orbital plane, the point in the orbit where the satellite passes the equator on its way up north-bound through the orbit is the ascending node, and the point where the satellite passes the equator on its way down south-bound is the descending node. Drawing a line through these two points on the equator defines the line of nodes.
Right ascension of ascending node (Ω) — The horizontal orientation of the ascending node of the ellipse (where the orbit passes upward through the reference plane) with respect to the I-axis of the reference frame. The rotation of the right ascension of the ascending node (RAAN) can be any number between 0° and 360°.

The remaining two elements are:

Argument of periapsis (ω) — The orientation of the ellipse in the orbital plane, as an angle measured from the ascending node to the periapsis in the range [0, 360).
True Anomaly (v) — The position of the orbiting body along the ellipse at a specific time. The position of the satellite on the path is measured counterclockwise from the periapsis and is called the true anomaly, ν, in the range [0, 360).

You can propagate the orbit using numerical methods such as numerical integration or Kepler's equation using these parameters. These methods allow you to calculate the position and velocity of the object at a given time. For information on each of these elements, see Orbital Elements.

The Aerospace Toolbox uses universal variables and Newton-Raphson iteration to propagate satellite orbits over time. This analytical algorithm is fast, but has limitations. Propagated orbits account only for central body spherical (point-mass) gravity. This formulation includes no other perturbations.

This propagation method is always performed in the ICRF inertial coordinate system with origin at the center of the central body. Given initial inertial position r₀ and velocity v₀ at time t₀, first find orbital energy, ξ, and the reciprocal of the semi-major axis, α:

$\begin{array}{l} ξ = \frac{v_{0}^{2}}{2} - \frac{μ}{r_{0}} \\ α = \frac{- 2 ξ}{μ}, \end{array}$

where μ is the standard gravitation parameter of the central body. Next, determine the orbit type from the sign of α.

α>0 => Circular or elliptical
α<0 => Hyperbolic
α≈0 => Parabolic

To initialize the Newton-Raphson iteration, select an initial guess for χ based on the orbit type:

Circular or elliptical orbit
$χ_{0} \approx \sqrt{μ} (Δ t) α,$
where Δt is the propagation step size (simulation time step). If Δt exceeds the orbital period $T = 2 π \sqrt{\frac{a^{3}}{μ}}$ , wrap Δt.
Parabolic orbit
$χ_{0} \approx \sqrt{p} 2 \cot (2 w),$
where:
$\begin{array}{l} \vec{h} = \vec{r_{0}} \times \vec{v_{0}} \\ p = \frac{h \cdot h}{μ} \\ \cot (2 s) = 3 \sqrt{\frac{μ}{p^{3}}} (Δ t) \\ \tan^{3} (w) = \tan (s) . \end{array}$
Hyperbolic orbit:
$χ_{0} \approx sign (Δ t) \sqrt{- \frac{1}{α}} \ln (\frac{- 2 μ α (Δ t)}{\vec{r_{0}} \cdot \vec{v_{0}} + sign (Δ t) \sqrt{- \frac{μ}{α}} (1 - r_{0} α)}) .$
Perform Newton-Raphson iteration while |x_n-x_n-1| > tolerance.
$\begin{array}{l} χ_{n + 1} = χ_{n} + \frac{\sqrt{μ} (Δ t) - χ_{n}^{3} c_{3} - \frac{\vec{r_{0}} \cdot \vec{v_{0}}}{\sqrt{μ}} χ_{n}^{2} c_{2} - r_{0} χ_{n} (1 - ψ c_{3})}{χ_{n}^{2} c_{2} + \frac{\vec{r_{0}} \cdot \vec{v_{0}}}{\sqrt{μ}} χ_{n} (1 - ψ c_{3}) + r_{0} (1 - ψ c_{2})} \\ χ_{n} \Leftarrow χ_{n + 1}, \end{array}$
where:
$ψ = χ_{n}^{2} α .$
(if ψ>0),
$\begin{array}{l} c_{2} = \frac{1 - \cos (\sqrt{ψ})}{ψ} \\ c_{3} = \frac{\sqrt{ψ} - \sin (\sqrt{ψ})}{\sqrt{ψ^{3}}} . \end{array}$
(if ψ<0),
$\begin{array}{l} c_{2} = \frac{1 - \cosh (\sqrt{- ψ})}{ψ} \\ c_{3} = \frac{\sinh (\sqrt{- ψ}) - \sqrt{- ψ}}{\sqrt{{(- ψ)}^{3}}} . \end{array}$
(if ψ≈0),
$\begin{array}{l} c_{2} = \frac{1}{2} \\ c_{3} = \frac{1}{6} . \end{array}$
Calculate universal variables $f$ , $\dot{f}$ , $g$ , and $\dot{g}$ .
$\begin{array}{l} f = 1 - \frac{χ_{n}^{2}}{r_{0}} c_{2} \\ \dot{f} = \frac{\sqrt{μ}}{r r_{0}} χ_{n} (ψ c_{3} - 1) \\ g = (Δ t) - \frac{χ_{n}^{3}}{\sqrt{μ}} c_{3} \\ \dot{g} = 1 - \frac{χ_{n}^{2}}{r} c_{2} . \end{array}$
Assemble position and velocity output vectors:
$\begin{array}{l} \vec{r_{icrf}} = f \vec{r_{0}} + g \vec{v_{0}} \\ \vec{v_{icrf}} = \dot{f} \vec{r_{0}} + \dot{g} \vec{v_{0}} . \end{array}$

Numerical (high precision)

This option uses the Simulink^® solver to integrate position and velocity from central body gravitational acceleration at each simulation timestep (Δt). The method for computing central body acceleration depends on the current setting for the Gravitational potential model parameter. You can also include custom acceleration components in the propagation algorithm using the A_icrf (applied acceleration) input port. For gravity models that include nonspherical acceleration terms, the block computes nonspherical gravity in a fixed-frame coordinate system (ITRF, in the case of Earth). Numerical integration, however, is always performed in the inertial ICRF coordinate system. Therefore, at each timestep:

The block transforms position and velocity states into the fixed-frame.
The block calculates nonspherical gravity in the fixed-frame.
The block transforms resulting acceleration into the inertial frame, where it is summed with the other acceleration terms and integrated.

$\begin{array}{l} \vec{a_{icrf}} = {\vec{a}}_{central body gravity} + {\vec{a}}_{drag} + {\vec{a}}_{thirdbody} + {\vec{a}}_{srp} + {\vec{a}}_{applied} \\ \vec{a_{icrf}} \underset{integrate}{\Rightarrow} \vec{r_{icrf}}, \vec{v_{icrf}} \end{array}$

where ${\vec{a}}_{applied}$

Central Body

Methods for computing central body acceleration include point-mass, oblate ellipsoid, or spherical harmonics.

Point-mass — This option treats the central body as a point-mass, including only the effects of spherical gravity using Newton's law of universal gravitation.
${\vec{a}}_{centralbodygravity} = - \frac{μ}{r^{2}} \frac{\vec{r_{icrf}}}{r},$
where μ is the standard gravitation parameter of the central body.
Oblate ellipsoid (J2) — In addition to spherical gravity, this option includes the perturbing effects of the second-degree, zonal harmonic gravity coefficient J₂, accounting for the oblateness of the central body. J₂ accounts for the vast majority of the central bodies gravitational departure from a perfect sphere.
${\overset{⇀}{a}}_{c e n t r a l b o d y g r a v i t y} = - \frac{μ}{r^{2}} \frac{\vec{r_{icrf}}}{r} + f i x e d 2 i n e r t i a l ({\overset{⇀}{a}}_{n o n s p h e r i c a l}),$
where:
$\begin{array}{l} {\vec{a}}_{nonspherical} = \\ {[\frac{1}{r} \frac{\partial}{\partial r} U - \frac{r_{ff}_{_{k}}}{r^{2} \sqrt{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}}} \frac{\partial}{\partial ϕ} U] r_{ff}_{_{i}}} i \\ + {[\frac{1}{r} \frac{\partial}{\partial r} U + \frac{r_{ff}_{_{k}}}{r^{2} \sqrt{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}}} \frac{\partial}{\partial ϕ} U] r_{ff}_{_{j}}} j \\ + {\frac{1}{r} (\frac{\partial}{\partial r} U) r_{k} + \frac{\sqrt{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}}}{r^{2}} \frac{\partial}{\partial ϕ} U} k, \end{array}$
given the partial derivatives in spherical coordinates:
$\begin{array}{l} \frac{\partial}{\partial r} U = \frac{3 μ}{r^{2}} (^{\frac{R_{cb}}{r}) 2} P_{2, 0} [\sin (ϕ)] J_{2} \\ \frac{\partial}{\partial ϕ} U = - \frac{μ}{r} (^{\frac{R_{cb}}{r}) 2} P_{2, 1} [\sin (ϕ)] J_{2} \end{array}$
where:
- ϕ and λ — Satellite geocentric latitude and longitude.
- P_2,0 and P_2,1 — Associated Legendre functions.
- μ — Standard gravitation parameter of the central body.
- R_cb — Central body equatorial radius.
The transformation fixed2inertial converts fixed-frame position, velocity, and acceleration into the ICRF coordinate system with origin at the center of the central body, accounting for centrifugal and coriolis acceleration. For more information about the fixed and inertial coordinate systems used for each central body, see Coordinate Systems for Navigation.
Spherical Harmonics — This option adds increased fidelity by including higher-order perturbation effects accounting for zonal, sectoral, and tesseral harmonics. For reference, the second-degree, zeroth order zonal harmonic J₂=-C₂₀. The Spherical Harmonics model accounts for harmonics up to max degree l=l_max, which varies by central body and geopotential model.
${\overset{⇀}{a}}_{c e n t r a l b o d y g r a v i t y} = - \frac{μ}{r^{2}} \frac{\vec{r_{icrf}}}{r} + f i x e d 2 i n e r t i a l ({\overset{⇀}{a}}_{n o n s p h e r i c a l}),$
where:
$\begin{array}{l} {\overset{⇀}{a}}_{n o n s p h e r i c a l} = \\ {[\frac{1}{r} \frac{\partial}{\partial r} U - \frac{r_{ff}_{_{k}}}{r^{2} \sqrt{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}}} \frac{\partial}{\partial ϕ} U] r_{ff}_{_{i}} - [\frac{1}{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}} \frac{\partial}{\partial λ} U] r_{ff}_{_{j}}} i \\ + {[\frac{1}{r} \frac{\partial}{\partial r} U + \frac{r_{ff}_{_{k}}}{r^{2} \sqrt{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}}} \frac{\partial}{\partial ϕ} U] r_{ff}_{_{j}} + [\frac{1}{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}} \frac{\partial}{\partial λ} U] r_{ff}_{_{i}}} j \\ + {\frac{1}{r} (\frac{\partial}{\partial r} U) r_{ff}_{_{k}} + \frac{\sqrt{r_{ff} {_{_{i}}}^{2} + r_{ff} {_{_{j}}}^{2}}}{r^{2}} \frac{\partial}{\partial ϕ} U} k, \end{array}$
given the following partial derivatives in spherical coordinates:
$\begin{array}{l} \frac{\partial}{\partial r} U = - \frac{μ}{r^{2}} {\sum^{}}_{l = 2}^{l_{\max}} {\sum^{}}_{m = 0}^{l} (^{\frac{R_{cb}}{r}) l} (l + 1) P_{l, m} [\sin (ϕ)] {C_{l, m} \cos (m λ) + S_{l, m} \sin (m λ)} \\ \frac{\partial}{\partial ϕ} U = \frac{μ}{r} {\sum^{}}_{l = 2}^{l_{\max}} {\sum^{}}_{m = 0}^{l} (^{\frac{R_{cb}}{r}) l} {P_{l, m + 1} [\sin (ϕ)] - (m) \tan (ϕ) P_{l, m} [\sin (ϕ)]} {C_{l, m} \cos (m λ) + S_{l, m} \sin (m λ)} \\ \frac{\partial}{\partial λ} U = \frac{μ}{r} {\sum^{}}_{l = 2}^{l_{\max}} {\sum^{}}_{m = 0}^{l} (^{\frac{R_{cb}}{r}) l} (m) P_{l, m} [\sin (ϕ)] {S_{l, m} \cos (m λ) - C_{l, m} \sin (m λ)}, \end{array}$
where:
- ϕ and λ — Satellite geocentric latitude and longitude.
- P_l,m — Associated Legendre functions.
- μ — Standard gravitation parameter of the central body.
- R_cb — Central body equatorial radius.
- C_l,m and S_l,m — Nonnormalized harmonic coefficients.
The transformation fixed2inertial converts fixed-frame position, velocity, and acceleration into the ICRF coordinate system with origin at the center of the central body, accounting for centrifugal and coriolis acceleration. For more information about the fixed and inertial coordinate systems used for each central body, see Coordinate Systems for Navigation.

Atmospheric Drag

Atmospheric drag is a force that opposes the motion of an object moving through a fluid, such as the atmosphere of the earth. Atmospheric drag is influenced by several factors, including the shape, size, and velocity of the object, and the properties of air. The primary mechanism behind atmospheric drag is the conversion of the kinetic energy of the object into heat through the process of air molecule collision and friction.

atmospheric drag of satellite

The Aerospace Toolbox uses this atmospheric drag equation:

$a_{d r a g} = - \frac{1}{2} ρ (\frac{C_{D} A}{m}) v^{2}_{r e l} \frac{{\vec{v}}_{r e l}}{v_{r e l}}$

where:

m — Spacecraft mass used by atmospheric drag calculation.
C_D — Coefficient of drag assuming that it is dimensionless.
ρ — Atmospheric density.
A — Area normal to v_rel, where
${\vec{v}}_{r e l} = {\vec{v}}_{s a t} + {\vec{v}}_{a t m o s}$
v_rel — Velocity relative to atmosphere.
${\vec{v}}_{r e l} = {\vec{v}}_{i c r f} - {\vec{ω}}_{\oplus} \times \vec{r_{i c r f}}$
where ${\vec{ω}}_{\oplus}$ is the central body angular velocity.

Third Body

In Aerospace Toolbox, third body refers to an additional celestial object that influences the motion of two primary bodies in a gravitational system.

This equation incorporates the third body contribution, enabling a more precise prediction of the motion of objects in space.

$a_{_{t h i r d_b o d y}} = μ_{t h i r d} (\frac{{\vec{r}}_{s a t, 3}}{r^{3}_{s a t, 3}} - \frac{\vec{r}}{r^{3}})$

where:

μ_third — Gravitational parameter of the third body.
${\vec{r}}_{s a t, 3}$ — Vector from the satellite to the third body.
$\vec{r}$ — Position of third body with regard to the central body, specified as a vector.

Solar Radiation Pressure

Solar Radiation Pressure (SRP) is the force produced by the impact of sunlight photons on the surface of an object in space. SRP is considered as a perturbing force that needs to be accounted for in orbit determination and prediction.

The Aerospace Toolbox uses this solar radiation pressure equation:

$a_{s r p} = - υ C_{r} \frac{A_{s}}{m} P_{s r p} {(\frac{A U}{r_{s a t, s u n}})}^{2} \frac{{\vec{r}}_{s a t, s u n}}{r_{s a t, s u n}}$

where:

m — Mass of the spacecraft
v — Eclipse shadow function
C_r — Spacecraft coefficient of reflectivity
A_s — Spacecraft solar radiation pressure area
P_srp — Solar radiation pressure at a distance of AU from the Sun
AU — Mean distance from the Sun to Earth (1AU)
| ${\vec{r}}_{_{s a t, s u n}}$ | — Vector from the satellite to the Sun origin

References

[1] Vallado, David, et al. “Revisiting Spacetrack Report #3.” AIAA/AAS Astrodynamics Specialist Conference and Exhibit, American Institute of Aeronautics and Astronautics, 2006. DOI.org (Crossref), https://doi.org/10.2514/6.2006-6753

Related Examples

Model, Visualize, and Analyze Satellite Scenario

More About

Satellite Scenario Key Concepts