Equations of Motion

The Equations of Motion eom library provides astrodynamics related functionality that can aid in the rapid development of analysis tools. The eomx application is included with the library. It is a command line program driven by a simple text based input file. From the command line, eomx is invoked:

eomx input_file.emx <iers_eop_file.txt>

The input_file.emx file defines the scenario and associated analysis tasks to execute. The file extension is arbitrary.

An optional IERS EOP (International Earth Rotation Service Earth Orientation Parameters) file can be specified via iers_eop_file.txt. When not included, IERS EOP values (polar motion, UT1-UTC, LOD, corrections to the reduction theory) are all set to zero. The .csv formatted finals.daily (IAU2000) or finals.all (IAU2000) files are accepted. Note values from these files are separated by semicolons, not commas. The first row is composed of column labels. The files are available at the top of the page from the version metadata links, not the unformatted latest version link.

Contents

Input File Overview
General Scenario Inputs
Date and Time
Duration
Units
Orbit Definition
Relative Orbit Definition
External Ephemerides
Ground Point Definition
Access Analysis to Ground Points
Commands

Example Input Case Files

Example Projects Built with eom via Makefile

State Vector Type

The state vector type indicates the type of coordinate system is used to specify state vector components.

CART

Cartesian coordinates composed of three position and three velocity elements.

KEP_T

Keplerian elements composed of semimajor axis, eccentricity, inclination, right ascension of the ascending node, argument of perigee, and true anomaly.

  AngleUnits  Degrees;
  # Default units of DU for semimajor axis
  #               a      e     i      o      w     v
  KEP_T  GCRF  4.1632  0.741  63.4  345.0  270.0  0.0;

Reference Frame

Not all settings taking a reference frame type accept all available reference frames.

GCRF

Cartesian geocentric celestial reference frame

ITRF

Cartesian international terrestrial reference frame

RTC

Cartesian radial, transverse, cross-track reference frame

Orbit Definition

Orbit definitions are named objects. They are defined by the Orbit keyword, the name of the object, epoch time, state vector, and propagator. General perturbation (GP) methods require no further inputs. Special Perturbation (SP) methods may include force model and integration method options. If not specified, the force model will default to two-body gravity with propagation via RK4 integration. Two examples will aid in understanding the input options to be described.

DistanceUnits  Kilometers;
TimeUnits Seconds;
Orbit  vinti    Vinti6   GD 2021 11 12 17 00 00.0
       CART  GCRF  -5552.0  -2563.0  3258.0   2.149  -7.539  -2.186;
Orbit  spj3     SP       GD 2021 11 12 17 00 00.0
       CART  GCRF  -5552.0  -2563.0  3258.0   2.149  -7.539  -2.186
       Propagator  RK4 Minutes 0.0;
       GravityModel  Jn 3;
       #Propagator  Adams4 Seconds 8.0;
       #GravityModel  Gravt 3 0;

The Orbit keyword is followed by a unique identifier, the orbit name. The orbit propagation model is next, followed by the epoch using a DateTime input option. The state vector type and reference frame follow. The final entry for GP propagation models is the state vector itself. The DistanceUnits and TimeUnits keywords affect these entries. SP methods can also include GravityModel include MoonGravity include SunGravity and Propagator (Numerical Integration) method.

Orbit Propagation Models

Kepler1

Kepler1 two-body GP model based off of Vinti's work and implemented by Gim J. Der and Herbert B. Reynolds, original available via AIAA.

Vinti6

Vinti6 GP model with full J2 and partial J3 zonal effects, based off Vinti's work and implemented by Gim J. Der and Herbert B. Reynolds, original available via AIAA. Unlike most analytic propagators, the Vinti model is initialized with state vectors, or those derived directly from osculating element sets (vs. mean element sets). For more information, see Orbital and Celestial Mechanics, Nino L. Bonavito, Gim J. Der, and John P. Vinti, AIAA, 1998.

VintiJ2

The Vinti6 model described above, but limited to only the J2 zonal spherical harmonic. A true J2 propagator (vs. the typical secular versions requiring mean element sets).

VintiMod

The Vinti6 model described above, but limited to only the J2 zonal spherical harmonic. A true J2 propagator (vs. the typical secular versions requiring mean element sets). This version differs from VintiJ2 in that it has been reworked, fully separating initialization from propagation, increasing computational efficiency.

SecJ2

GENPL option implementing standard secular J2 effects, best initialized with mean element sets.

OscJ2

GENPL option implementing a highly optimized version of the Vinti propagator.

SP

Special perturbations methods. The equations of motion are resolved through the use of numerical integration.

Gravity Model

SP methods can select the gravity model to be used. The EGM 2008 gravity model is employed for all implementations. GravityModel Jn 0 is the default if not specified.

GravityModel Jn d

A simple zonal gravity model supporting up to degree (d) 4

GravityModel Standard d o

Standard eom rectangular gravity model supporting degree and order d o such that order .LE. degree. The implementation is based primarily off equations from Vallado's Fundamentals of Astrodynamics and Applications, 3rd Ed.

GravityModel Gravt d o

GENPL option supporting degree and order d o

Moon Gravity Model

SP methods can select the model to be used for gravitational perturbations due to the moon. The default is to not incorporate lunar gravitational effects.

MoonGravity Meeus

Moon gravitational perturbations with lunar ephemeris resolved through the method from Chapter 47 of Meeus' "Astronomical Algorithms", 2nd Ed.

Sun Gravity Model

SP methods can select the model to be used for gravitational perturbations due to the sun. The default is to not incorporate solar gravitational effects.

SunGravity Meeus

Sun gravitational perturbations with solar ephemeris resolved through the method from Chapter 25 of Meeus' "Astronomical Algorithms", 2nd Ed.

Numerical Integrator

SP methods can select the numerical integration method, used to solve the equations of motion. Propagator RK4 is the default if not specified.

Propagator RK4 Duration

Standard Runge-Kutta integration method. Robust, but not efficient. The Duration is the integration step size to use. A Duration equivalent to zero will result in the use of a default integration step size.

Propagator Adams4 Duration

Fixed step size Adams-Bashforth predictor with Adams-Moulton corrector, primed via RK4. The Duration is the integration step size to use. A Duration equivalent to zero will result in the use of a default integration step size.

Propagator GJ Duration

GENPL option, highly efficient Gauss-Jackson 2nd order integrator with on-the-fly step size adjustment to minimize error while maximizing speed. The Duration entry is required by the parser, but the setting is currently ignored. Instead, a default initial integration step size is used and adjusted during propagation.

Propagator GJs Duration

GENPL option, GJ propagator employing time regularization. Better suited for elliptical orbits. As with the GJ propagation, Duration is ignored.

Relative Orbit Definition

Allows for a relative orbit to be specified based on a RTC relative bounding box and transverse distance offset. The method is based on Spacecraft Relative Orbit Geometry Description Through Orbit Element Differences by Hanspeter Schaub. This method of defining a bounding box and offset automatically guarantees the energy matching constraint. The resulting orbit definitions are propagated in inertial space, not a linearized reference frame intended for control systems only. This method is not limited to circular orbits.

The syntax is:

RelativeOrbit orbit_name template_orbit_name RTC dR dT dC deltaT

The values dR dT dC define the desired bounding box limits of the relative orbit. Some values may be exceeded to meet minimum values. The deltaT parameter is an along-track offset between the orbit centers.

External Ephemerides

Externally generated SP3c is the only format currently supported for ingestion. Examples of loading SP3c formatted ephemeris from the file nsgf.orb.lageos2.210626.v70.sp3 with orbit names nsgf_lageos2h and nsgf_lageos2t are:

EphemerisFile nsgf_lageos2h  SP3c Hermite nsgf.orb.lageos2.210626.v70.sp3;
EphemerisFile nsgf_lageos2t  SP3c Chebyshev nsgf.orb.lageos2.210626.v70.sp3;

Hermite interpolation is performed for the first entry. Each interpolator is fit to two state vectors (position and velocity) with acceleration generated by a simple J4 gravity model. Hermite interpolation assumes velocity is the derivative of position, and should therefore be suitable for externally generated precision ephemerides. This assumption would not necessarily be valid for ephemererides generated via analytic (GP) propagation methods. The assumption that a J4 gravity model is sufficient to improve the fit must also be valid (typically a very good assumption). The resulting 5th order polynomial explicitly maintains continuity through acceleration and smoothness through velocity. Continuity and smoothness are implicitly maintained to higher degrees than the 2nd derivative. Chebyshev interpolation currently fits position and velocity separately, eliminating the requirement that velocity be the derivative of position. This interpolation method currently fits 9 points to an 8th order polynomial, ensuring endpoints match between polynomials. More refinement to the implementation is expected, most likely through the fitting of more points with endpoint matching constraints.

Ground Point Definition

Ground points can be defined given geodetic coordinates latitude, longitude, and altitude (w.r.t. the ellipsoid). Units are based on default or preceding Sticky entries.

GroundPoint name geodetic_latitude longitude altitude

Access Analysis to Ground Points

Access analysis from a ground point to an orbit can be defined with the GroundPoint keyword:

Access GroundPointAccess AccessModel orbit_name ground_point_name

Two access models are supported. The Standard model uses a variable step size based on the rate of change in true anomaly to locate the next access window. Once the access window is located, a bisection method is used to refine the precise rise and set times. This method is robust and fast. A significantly slower Debug option is also available that uses a small fixed step size to search for access windows. By default, access is generated based on line if sight visibility from the ground point to the satellite. This occurs when the elevation of the satellite w.r.t. the ground point, relative to the plane tangent to the ellipsoid, is zero. Currently, only the minimum elevation is configurable through an additional line:

MinimumElevation d

The following illustrates creating a ground point location named gp01 where positive latitude implies north and positive longitude implies east. Access to this location by a satellite named leosat is to be performed with a minimum ground elevation constraint of -5 degrees.

AngleUnits  Degrees;
GroundPoint gp01 LLA  43.0  30.0  0.0;
Access  GroundPointAccess  Standard  leosat gp01
                           MinimumElevation  -5.0;

Commands

Commands typially operate on named objects. Most are affected by the sticky OutputRate keyword. Examples:

Command PrintEphemeris gj_40x40 GCRF gj_40x40_i.e;
Command PrintEphemeris gj_40x40 ITRF gj_40x40_f.e;
OutputRate Minutes 1;
Command PrintRange vinti kep vinti_rng_vk;
Command PrintRange vinti vintij2 vinti_rng_vvj2;
Command PrintRTC vinti spj3 vinti_rtc_vsp;

OutputRate Duration

Sticky setting specifying the rate at which an output should be generated.

Command PrintEphemeris orbit_name frame output_file_name

Writes STK .e formatted position and velocity ephemeris for the specified satellite name. The output rate is determined internally. The filename suffix is user specified (does not default to .e). Output units are meters and meters/second.

Command PrintRange orbit1_name orbit2_name output_file_prefix

Writes a Matlab/Octave compatible function to the specified filename with the .m suffix added. When run in Matlab, the function plots the range vs. time over the duration of the scenario. The figure handle and data are returned for use. Output units are specified by DistanceUnits.

Command PrintRangeSpectrum orbit1_name orbit2_name output_file_prefix

Writes a Matlab/Octave compatible function to the specified filename with the .m suffix added. When run in Matlab, the function plots the range vs. time over the duration of the scenario. The figure handle and data are returned for use. Output units are specified by DistanceUnits. In addition, the spectrum of the error is plotted where frequency is revs/period w.r.t. orbit1. This functionality is intended for comparing differences between ephemerides of the same orbit definition. For example, to compare the error between SLR based and numerically propagated ephemerides. Or, say, the difference between two analytic propagators.

Command PrintRTC orbit1_name orbit2_name output_file_prefix

Writes a Matlab/Octave compatible function to the specified filename with the .m suffix added. When run in Matlab, the function plots the relative position of the second orbit w.r.t. the first, vs. time over the duration of the scenario. The 3D plot is in the traditional radial, transverse (in-track), cross-track reference frame used for relative orbital dynamics. The figure handle and data are returned for use. Output units are specified by DistanceUnits.