https://github.com/armavica/sgp4
A Rust implementation of the SGP4 satellite propagation algorithm
Science Score: 13.0%
This score indicates how likely this project is to be science-related based on various indicators:
-
○CITATION.cff file
-
○codemeta.json file
-
○.zenodo.json file
-
✓DOI references
Found 6 DOI reference(s) in README -
○Academic publication links
-
○Academic email domains
-
○Institutional organization owner
-
○JOSS paper metadata
-
○Scientific vocabulary similarity
Low similarity (8.4%) to scientific vocabulary
Last synced: 10 months ago
·
JSON representation
Repository
A Rust implementation of the SGP4 satellite propagation algorithm
Basic Info
- Host: GitHub
- Owner: Armavica
- License: mit
- Default Branch: master
- Homepage: https://crates.io/crates/sgp4
- Size: 339 KB
Statistics
- Stars: 0
- Watchers: 1
- Forks: 0
- Open Issues: 0
- Releases: 0
Fork of neuromorphicsystems/sgp4
Created about 5 years ago
· Last pushed about 5 years ago
https://github.com/Armavica/sgp4/blob/master/
# SGP4
The SGP4 algorithm, ported to Rust from the reference Celestrak implementation [[1]](#1).
The code was entirely refactored to leverage Rust's algebraic data types and highlight the relationship between the implementation and the reference mathematical equations [[2]](#2).
The numerical predictions are almost identical to those of the Celestrak implementation. The observed differences (less than 2 10 km for the position and 10 km.s for the velocity three and a half years after epoch) are well below the accuracy of the algorithm.
We drew inspiration from the incomplete https://github.com/natronics/rust-sgp4 to write mathematical expressions using UTF-8 characters.
## Cargo
```
[dependencies]
sgp4 = "0.5"
```
## Documentation
The code documentation is hosted at [https://docs.rs/sgp4/0.5.0/sgp4/](https://docs.rs/sgp4/0.5.0/sgp4/).
Examples can be found in this repository's *examples* directory:
- *examples/celestrak.rs* retrieves the most recent "stations" OMMs from Celestrak and propagates them
- *examples/omm.rs* parses and propagates a JSON-encoded OMM
- *examples/space-track.rs* retrieves the 20 most recent launches OMMs from Space-Track and propagates them
- *examples/tle.rs* parses and propagates a TLE
- *examples/tle_afspc.rs* parses and propagates a TLE using the AFSPC compatibility mode
- *examples/advanced.rs* leverages the advanced API to (marginally) accelerate the propagation of deep space resonant satellites
To run an example (here *examples/celestrak.rs*), use:
```sh
cargo run --example celestrak
```
To run the Space-Track example, you must first assign your Space-Track.org credentials to the fields `identity` and `password` (see lines 3 and 4 in *examples/space-track.rs*).
## Benchmark
The benchmark code is available at https://github.com/neuromorphicsystems/sgp4-benchmark. It compares two SGP4 implementations in different configurations:
- `cpp`: the Celestrak implementation [[1]](#1) in improved mode
- `cpp-afspc`: the Celestrak implementation [[1]](#1) in AFSPC compatibility mode
- `cpp-fastmath`: the Celestrak implementation [[1]](#1) in improved mode with the `fast-math` compiler flag
- `cpp-afspc-fastmath`: the Celestrak implementation [[1]](#1) in AFSPC compatibility mode with the `fast-math` compiler flag
- `rust`: our Rust implementation in default mode
- `rust-afspc`: our Rust implementation in AFSPC compatibility mode
This benchmark must not be confused with the code in this repository's *bench* directory. The latter considers only a small subset of the Celestrak catalogue (the tests recommended in [[1]](#1)) and does not measure the original C++ implementation.
The present results were obtained using a machine with the following configuration:
- __CPU__ - Intel Core i7-8700 @ 3.20GHz
- __RAM__ - Kingston DDR4 @ 2.667 GHz
- __OS__ - Ubuntu 16.04
- __Compilers__ - Rust 1.44.1 and gcc 9.3.0
Accuracy measures the maximum propagation error of each implementation with respect to the reference implementation (`cpp-afspc`) over the full Celestrak catalogue (1 minute timestep over 24 hours).
| implementation | maximum position error | maximum speed error |
|----------------------|------------------------|---------------------|
| `cpp-afspc` | (reference) | (reference) |
| `cpp` | 1.05 km | 1.30 10 km.s |
| `cpp-fastmath` | 1.05 km | 1.30 10 km.s |
| `cpp-afspc-fastmath` | 4.21 10 km | 7.51 10 km.s |
| `rust` | 1.05 km | 1.30 10 km.s |
| `rust-afspc` | 4.19 10 km | 7.46 10 km.s |
The Rust and C++ fast-math errors have the same order of magnitude. In both cases, they can be attributed to mathematically identical expressions implemented with different floating-point operations.
Speed measures the time it takes to propagate every satellite in the Celestrak catalogue (1 minute timestep over 24 hours) using a single thread. 100 values are sampled per implementation.
| implementation | minimum | Q1 | median | Q3 | maximum | relative difference |
|----------------------|---------|--------|--------|--------|---------|---------------------|
| `cpp-afspc` | 8.95 s | 9.02 s | 9.03 s | 9.06 s | 9.18 s | (reference) |
| `cpp` | 8.95 s | 9.01 s | 9.04 s | 9.06 s | 9.25 s | + 0 % |
| `cpp-fastmath` | 7.67 s | 7.74 s | 7.77 s | 7.79 s | 7.90 s | - 14 % |
| `cpp-afspc-fastmath` | 7.70 s | 7.74 s | 7.76 s | 7.79 s | 7.86 s | - 14 % |
| `rust` | 8.36 s | 8.41 s | 8.43 s | 8.45 s | 8.53 s | - 7 % |
| `rust-afspc` | 8.36 s | 8.41 s | 8.43 s | 8.46 s | 8.59 s | - 7 % |
Rust fast-math support is a work in progress (see https://github.com/rust-lang/rust/issues/21690). Similarly to C++, it should have a very small impact on accuracy while providing a substantial speed gain.
## Variables and mathematical expressions
### Variables
Each variable is used to store the result of one and only one expression. Most variables are immutable, with the exception of the variable `(E + )` used to solve Kepler's equation and the state variables `t`, `n` and `` used to integrate the resonance effects of Earth gravity.
The following tables list the variables used in the code and their associated mathematical symbol. Where possible, we used symbols from [[2]](#2). Partial expressions without a name in [[2]](#2) follow the convention `k, n ` if they are shared between initialization and propagation, and `p, n ` if they are local to initialization or propagation.
1. [Initialization variables](#initialization-variables)
2. [Propagation variables](#propagation-variables)
3. [Third-body initialization variables](#third-body-initialization-variables)
4. [Third-body propagation variables](#third-body-propagation-variables)
---
#### Initialization variables
The following variables depend solely on epoch elements.
| variable | symbol | description |
|:----------------------------------------|:---------------|:------------|
| `Elements::datetime.year()` | `y` | Gregorian calendar year |
| `Elements::datetime.month()` | `m` | Gregorian calendar month in the range `[1, 12]` |
| `Elements::datetime.day()` | `d` | Gregorian calendar day in the range `[1, 31]` |
| `Elements::datetime.hour()` | `h` | Hours since midnight in the range `[0, 23]` |
| `Elements::datetime.minute()` | `min` | Minutes since the hour in the range `[0, 59]` |
| `Elements::datetime.second()` | `s` | Seconds since the minute in the range `[0, 59]` |
| `Elements::datetime.nanosecond()` | `ns` | Nanoseconds since the second in the range `[0, 10[` |
| `epoch` | `y` | Julian years since UTC 1 January 2000 12h00 (J2000) |
| `d1900` | `d` | Julian days since UTC 1 January 1900 12h00 (J1900) |
| `d1970` | `d` | Julian days since UTC 1 January 1970 12h00 (J1970) |
| `c2000` | `c` | Julian centuries since UTC 1 January 2000 12h00 (J2000) |
| `geopotential.ae` | `a` | equatorial radius of the earth in km |
| `geopotential.ke` | `k` | square root of the earth's gravitational parameter in earth radii min |
| `geopotential.j2` | `J` | un-normalised second zonal harmonic |
| `geopotential.j3` | `J` | un-normalised third zonal harmonic |
| `geopotential.j4` | `J` | un-normalised fourth zonal harmonic |
| `kozai_mean_motion` | `n` | mean number of orbits per day (Kozai convention) at epoch in rad.min |
| `a1` | `a` | semi-major axis at epoch (Kozai convention) |
| `p0` | `p` | partial expression of `` and `` |
| `d1` | `` | used in the Kozai to Brouwer conversion |
| `d0` | `` | used in the Kozai to Brouwer conversion |
| `B*` | `B*` | radiation pressure coefficient in earth radii |
| `orbit_0.inclination` | `I` | angle between the equator and the orbit plane at epoch in rad |
| `orbit_0.right_ascension` | `` | angle between vernal equinox and the point where the orbit crosses the equatorial plane at epoch in rad |
| `orbit_0.eccentricity` | `e` | shape of the orbit at epoch |
| `orbit_0.argument_of_perigee` | `` | angle between the ascending node and the orbit's point of closest approach to the earth at epoch in rad |
| `orbit_0.mean_anomaly` | `M` | angle of the satellite location measured from perigee at epoch in rad |
| `orbit_0.mean_motion` | `n"` | mean number of orbits per day (Brouwer convention) at epoch in rad.min |
| `p1` | `p` | cosine of the inclination at epoch used in multiple expressions during initialization (`` in [[2]](#2), renamed to avoid confusion with the sidereal time) |
| `p2` | `p` | partial expression of multiple initialization expressions |
| `a0` | `a"` | semi-major axis at epoch (Brouwer convention) |
| `p3` | `p` | perigee in earth radii |
| `p4` | `p` | height of perigee in km |
| `p5` | `p` | partial expression of `s` |
| `s` | `s` | altitude parameter of the atmospheric drag expression |
| `p6` | `p` | partial expression of the atmospheric drag |
| `xi` | `` | partial expression of multiple initialization expressions |
| `p7` | `p` | partial expression of multiple initialization expressions |
| `eta` | `` | partial expression of multiple initialization expressions and of the argument of perigee and mean anomaly in eccentric high altitude near earth propagation |
| `p8` | `p` | partial expression of multiple initialization expressions |
| `p9` | `p` | partial expression of multiple initialization expressions |
| `c1` | `C` | partial expression of multiple initializationa and propagation expressions |
| `p10` | `p` | partial expression of multiple initialization expressions |
| `b0` | `` | partial expression of multiple initialization expressions |
| `p11` | `p` | partial expression of multiple initialization expressions |
| `p12` | `p` | partial expression of multiple initialization expressions |
| `p13` | `p` | partial expression of multiple initialization expressions |
| `p14` | `p` | partial expression of multiple initialization expressions |
| `p15` | `p` | partial expression of multiple initialization expressions |
| `k14` | `k` | first order coefficient of the argument of perigee before adding solar and lunar perturbations |
| `c4` | `C` | partial expression of multiple initializationa and propagation expressions (differs from the `C` expression in [[2]](#2) by a factor B*) |
| `right_ascension_dot` | `` | first order coefficient of the right ascension |
| `argument_of_perigee_dot` | `` | first order coefficient of the argument of perigee |
| `mean_anomaly_dot` | `` | first order coefficient of the mean anomaly |
| `k0` | `k` | second order coefficient of the right ascension before adding perturbations |
| `k1` | `k` | partial expression of the second order coefficient of the mean anomaly |
| `k2` | `k` | partial expression of `a` in near earth propagation |
| `k3` | `k` | partial expression of `r`, `r` and `rf` in near earth propagation |
| `k4` | `k` | partial expression of `u` in near earth propagation |
| `k5` | `k` | partial expression of the initial Kepler variable `p` in near earth propagation |
| `k6` | `k` | partial expression of multiple initialization expressions and of `r` and `rf` in near earth propagation |
| `d2` | `D` | partial expression of multiple near earth initialization expressions and of the semi-major axis in near earth propagation |
| `p16` | `p` | partial expression of multiple near earth initialization expressions |
| `d3` | `D` | partial expression of multiple near earth initialization expressions and of the semi-major axis in near earth propagation |
| `d4` | `D` | partial expression of multiple near earth initialization expressions and of the semi-major axis in near earth propagation |
| `c5` | `C` | partial expression of multiple initializationa and propagation expressions (differs from the `C` expression in [[2]](#2) by a factor B*)
| `k7` | `k` | sine of the mean anomaly at epoch |
| `k8` | `k` | partial expression of the mean anomaly third order coefficient in high altitude near earth propagation |
| `k9` | `k` | partial expression of the mean anomaly fourth order coefficient in high altitude near earth propagation |
| `k10` | `k` | partial expression of the mean anomaly fifth order coefficient in high altitude near earth propagation |
| `k11` | `k` | partial expression of the argument of perigee and mean anomaly in eccentric high altitude near earth propagation |
| `k12` | `k` | partial expression of the argument of perigee and mean anomaly in eccentric high altitude near earth propagation |
| `k13` | `k` | partial expression of the argument of perigee and mean anomaly in eccentric high altitude near earth propagation |
| `lunar_right_ascension_epsilon` | `` | lunar right ascension of the ascending node |
| `lunar_right_ascension_sine` | `sin ` | sine of the lunar right ascension of the ascending node referred to the equator |
| `lunar_right_ascension_cosine` | `cos ` | cosine of the lunar right ascension of the ascending node referred to the equator |
| `lunar_argument_of_perigee` | `` | lunar argument of perigee |
| `sidereal_time_0` | `` | Greenwich sidereal time at epoch |
| `lambda_0` | `` | Earth gravity resonance variable at epoch |
| `lambda_dot_0` | `` | time derivative of the Earth gravity resonance variable at epoch |
| `p17` | `p` | partial expression of ``, `` and `` |
| `dr1` | `` | first Earth gravity resonance coefficient for geosynchronous orbits (`` in [[2]](#2), renamed to avoid confusion with `` used in the Kozai to Brouwer conversion) |
| `dr2` | `` | second Earth gravity resonance coefficient for geosynchronous orbits (`` in [[2]](#2), renamed to match ``) |
| `dr3` | `` | third Earth gravity resonance coefficient for geosynchronous orbits (`` in [[2]](#2), renamed to match ``) |
| `p18` | `p` | partial expression of `D` and `D` |
| `p19` | `p` | partial expression of `D` and `D` |
| `p20` | `p` | partial expression of `D` and `D` |
| `p21` | `p` | partial expression of `D`, `D`, `D` and `D` |
| `f220` | `F` | partial expression of `D` and `D` |
| `g211` | `G` | partial expression of `D` |
| `g310` | `G` | partial expression of `D` |
| `g322` | `G` | partial expression of `D` |
| `g410` | `G` | partial expession of `D` |
| `g422` | `G` | partial expession of `D` |
| `g520` | `G` | partial expression of `D` |
| `g532` | `G` | partial expression of `D` |
| `g521` | `G` | partial expression of `D` |
| `g533` | `G` | partial expression of `D` |
| `d2201` | `D` | gravity resonance coefficient for Molniya orbits (the `D` expression in [[2]](#2) is missing a factor `l - 2p + k` from the original equation in [[4]](#4) with `k = -1` instead of `1`) |
| `d2211` | `D` | gravity resonance coefficient for Molniya orbits (the `D` expression in [[2]](#2) is missing a factor `l - 2p + k` from the original equation in [[4]](#4)) |
| `d3210` | `D` | see `D` |
| `d3222` | `D` | see `D` |
| `d4410` | `D` | see `D` |
| `d4422` | `D` | see `D` |
| `d5220` | `D` | see `D` |
| `d5232` | `D` | see `D` |
| `d5421` | `D` | see `D` |
| `d5433` | `D` | see `D` |
#### Propagation variables
The following expressions depend on the propagation time `t`.
| variable | symbol | description |
|:----------------------------------------|:---------------|:------------|
| `t` | `t` | minutes elapsed since epoch (can be negative) |
| `p22` | `p` | right ascension of the ascending node with neither Earth gravity resonance nor Sun and Moon contributions |
| `p23` | `p` | argument of perigee with neither high altitude drag effects, Earth gravity resonance nor Sun and Moon contributions |
| `orbit.inclination` | `I` | inclination at epoch plus `t` without the short-period effects of Earth gravity |
| `orbit.right_ascension` | `` | right ascension of the ascending node at epoch plus `t` without the short-period effects of Earth gravity |
| `orbit.eccentricity` | `e` | eccentricity at epoch plus `t` without the short-period effects of Earth gravity |
| `orbit.argument_of_perigee` | `` | argument of perigee at epoch plus `t` without the short-period effects of Earth gravity |
| `orbit.mean_anomaly` | `M` | mean anomaly at epoch plus `t` without the short-period effects of Earth gravity |
| `orbit.mean_motion` | `n` | mean motion at epoch plus `t` without the short-period effects of Earth gravity |
| `a` | `a` | semi-major axis |
| `p32` | `p` | partial expression of `a` |
| `p33` | `p` | partial expression of `r`, `r` and `rf` |
| `p34` | `p` | partial expression of `u` |
| `p35` | `p` | partial expression of the initial Kepler variable `p` |
| `p36` | `p` | partial expression of `r` and `rf` |
| `p37` | `p` | partial expression of `a` and the initial Kepler variable `p` |
| `axn` | `a` | normalized linear eccentricity projected on the line of nodes |
| `ayn` | `a` | normalized linear eccentricity projected on the normal to the line of nodes |
| `p38` | `p` | initial Kepler variable (`U` in [[2]](#2), renamed to avoid confusion with the true anomaly plus argument of perigee `u`) |
| `ew` | `(E + )` | Kepler variable used in an iterative process to estimate the eccentric anomaly `E` |
| `delta ` | `(E + )` | correction to the Kepler variable at iteration `i` |
| `p39` | `p` | eccentricity at epoch plus `t` |
| `pl` | `p` | semi-latus rectum |
| `p40` | `p` | normalized linear eccentricity projected on the semi-minor axis |
| `r` | `r` | radius (distance to the focus) without the short-period effects of Earth gravity |
| `r_dot` | `r` | radius time derivative without the short-period effects of Earth gravity |
| `b` | `` | semi-minor axis over semi-major axis |
| `p41` | `p` | partial expression of `p` and `p` |
| `p42` | `p` | sine of `u` |
| `p43` | `p` | cosine of `u` |
| `u` | `u` | true anomaly plus argument of perigee without the short-period effects of Earth gravity |
| `p44` | `p` | `sin(2 u)`, partial expression of `u`, `` and `r` |
| `p45` | `p` | `cos(2 u)`, partial expression of `r`, `I` and `rf` |
| `p46` | `p` | partial expression of `r`, `u`, `I` and `` |
| `rk` | `r` | radius (distance to the focus) |
| `uk` | `u` | true anomaly plus argument of perigee |
| `inclination_k` | `I` | inclination at epoch plus `t`
| `right_ascension_k` | `` | right ascension at epoch plus `t`
| `rk_dot` | `r` | radius time derivative (@DEV orthogonal speed) |
| `rfk_dot` | `rf` | radius times the true anomaly derivative (@DEV non-orthogonal speed) |
| `u0` | `u` | x component of the position unit vector |
| `u1` | `u` | y component of the position unit vector |
| `u2` | `u` | z component of the position unit vector |
| `prediction.position[0]` | `r` | x component of the position vector in km (True Equator, Mean Equinox (TEME) of epoch reference frame) |
| `prediction.position[1]` | `r` | y component of the position vector in km (True Equator, Mean Equinox (TEME) of epoch reference frame) |
| `prediction.position[2]` | `r` | z component of the position vector in km (True Equator, Mean Equinox (TEME) of epoch reference frame) |
| `prediction.velocity[0]` | `r` | x component of the velocity vector in km.s (True Equator, Mean Equinox (TEME) of epoch reference frame) |
| `prediction.velocity[1]` | `r` | y component of the velocity vector in km.s (True Equator, Mean Equinox (TEME) of epoch reference frame) |
| `prediction.velocity[2]` | `r` | z component of the velocity vector in km.s (True Equator, Mean Equinox (TEME) of epoch reference frame) |
| `p24` | `p` | mean anomaly without drag contributions in near earth propagation |
| `p25` | `p` | partial expression of `` and `M` in near earth propagation |
| `p26` | `p` | mean anomaly with elliptic correction and without drag contributions in near earth propagation |
| `p27` | `p` | non-clamped eccentricity in near earth propagation |
| `p28` | `p` | semi-major axis with resonance correction in deep space propagation |
| `p29` | `p` | mean anomaly with resonance correction in deep space propagation |
| `p31` | `p` | non-clamped eccentricity in deep space propagation |
| `sidereal_time` | `` | sidereal time at epoch plus `t` |
| `delta_t` | `t` | time step used in the integration of resonance effects of Earth gravity in min (either `720` or `-720`) |
| `lambda_dot` | `` | resonance effects of Earth gravity variable's time derivative at epoch plus `i t` |
| `ni_dot` | `n` | mean motion time derivative at epoch plus `i t` |
| `ni_ddot` | `n` | mean motion second time derivative at epoch plus `i t` |
| `ResonanceState::t` | `t` | resonance effects of Earth gravity integrator time (`i t`) |
| `ResonanceState::mean_motion` | `n` | mean motion time derivative at epoch plus `t i` |
| `ResonanceState::lambda` | `` | resonance effects of Earth gravity variable at epoch plus `i t` |
| `p30` | `p` | non-normalised `` in Lyddane deep space propagation |
#### Third-body initialization variables
The contribution of the Sun and the Moon to the orbital elements are calculated with a unique set of expressions. *src/third_body.rs* provides a generic implementation of these expressions. Variables specific to the third body (either the Sun or the Moon) are annotated with `x`. In every other file, these variables are annotated with `s` if they correspond to solar perturbations, and `l` if they correspond to lunar perturbations.
The `a`, `X`, `Z` (`n `), `F` and `F` variables correspond to the `a`, `X`, `Z`, `F` and `F` variables in [[2]](#2). The added `x` highlights the dependence on the perturbing third body.
The following variables depend solely on epoch elements.
| variable | symbol | description |
|:----------------------------------------|:---------------|:------------|
| `third_body_inclination_sine` | `sin I` | sine of the inclination of the Sun (`sin I`) or the Moon (`sin I`) |
| `third_body_inclination_cosine` | `cos I` | cosine of the inclination of the Sun (`cos I`) or the Moon (`cos I`) |
| `delta_right_ascension_sine` | `sin( - )` | sine of the difference between the right ascension of the ascending node of the satellite at epoch and the Sun's (`sin( - )`) or the Moon's (`sin( - )`) |
| `delta_right_ascension_cosine` | `cos( - )` | cosine of the difference between the right ascension of the ascending node of the satellite at epoch and the Sun's (`cos( - )`) or the Moon's (`cos( - )`) |
| `third_body_argument_of_perigee_sine` | `sin ` | sine of the argument of perigee of the Sun (`sin `) or the Moon (`sin `) |
| `third_body_argument_of_perigee_cosine` | `cos ` | cosine of the argument of perigee of the Sun (`sin `) or the Moon (`cos `) |
| `third_body_mean_anomaly_0` | `M` | mean anomaly at epoch of the Sun (`M`) or the Moon (`M`) |
| `ax1` | `a` | partial expression of multiple `X` and `Z` expressions |
| `ax3` | `a` | partial expression of multiple `X` and `Z` expressions |
| `ax7` | `a` | partial expression of multiple `a` and `a` |
| `ax8` | `a` | partial expression of multiple `a` and `a` |
| `ax9` | `a` | partial expression of multiple `a` and `a` |
| `ax10` | `a` | partial expression of multiple `a` and `a` |
| `ax2` | `a` | partial expression of multiple `X` and `Z` expressions |
| `ax4` | `a` | partial expression of multiple `X` and `Z` expressions |
| `ax5` | `a` | partial expression of multiple `X` and `Z` expressions |
| `ax6` | `a` | partial expression of multiple `X` and `Z` expressions |
| `xx1` | `X` | partial expression of multiple `Z` expressions, `k`, `k` and `` |
| `xx2` | `X` | partial expression of multiple `Z` expressions, `k`, `k` and `` |
| `xx3` | `X` | partial expression of multiple `Z` expressions, `k`, `k` and `` |
| `xx4` | `X` | partial expression of multiple `Z` expressions, `k`, `k` and `` |
| `xx5` | `X` | partial expression of multiple `Z` expressions |
| `xx6` | `X` | partial expression of multiple `Z` expressions |
| `xx7` | `X` | partial expression of multiple `Z` expressions |
| `xx8` | `X` | partial expression of multiple `Z` expressions |
| `zx31` | `Z` | partial expression of `Z`, `k` and `` |
| `zx32` | `Z` | partial expression of `Z`, `k` and `` |
| `zx33` | `Z` | partial expression of `Z`, `k` and `` |
| `zx11` | `Z` | partial expression of `k` and `I` |
| `zx13` | `Z` | partial expression of `k` and `I` |
| `zx21` | `Z` | partial expression of `k` and `` |
| `zx23` | `Z` | partial expression of `k` and `` |
| `zx1` | `Z` | partial expression of `k` and `` |
| `zx3` | `Z` | partial expression of `k` and `` |
| `px0` | `p` | partial expression of multiple `k` expressions and `` |
| `px1` | `p` | partial expression of multiple `k` expressions and `I` |
| `px2` | `p` | partial expression of multiple `k` expressions and `` |
| `px3` | `p` | partial expression of multiple `k` expressions and `` |
| `kx0` | `k` | `F` coefficient of `e` |
| `kx1` | `k` | `F` coefficient of `e` |
| `kx2` | `k` | `F` coefficient of `I` |
| `kx3` | `k` | `F` coefficient of `I` |
| `kx4` | `k` | `F` coefficient of `M` |
| `kx5` | `k` | `F` coefficient of `M` |
| `kx6` | `k` | `sin f` coefficient of `M` |
| `kx7` | `k` | `F` coefficient of `p` |
| `kx8` | `k` | `F` coefficient of `p` |
| `kx9` | `k` | `sin f` coefficient of `p` |
| `kx10` | `k` | `F` coefficient of `p` |
| `kx11` | `k` | `F` coefficient of `p` |
| `third_body_dots.inclination` | `I` | secular contribution of the Sun (`I`) or the Moon (`I`) to the inclination |
| `third_body_right_ascension_dot` | `` | secular contribution of the Sun (``) or the Moon (``) to the right ascension of the ascending node |
| `third_body_dots.eccentricity` | `` | secular contribution of the Sun (``) or the Moon (``) to the eccentricity |
| `third_body_dots.agument_of_perigee` | `` | secular contribution of the Sun (``) or the Moon (``) to the argument of perigee |
| `third_body_dots.mean_anomaly` | `` | secular contribution of the Sun (``) or the Moon (``) to the mean anomaly |
#### Third-body propagation variables
The following variables depend on the propagation time `t`.
| variable | symbol | description |
|:----------------------------------------|:---------------|:------------|
| `third_body_mean_anomaly` | `M` | mean anomaly of the Sun (`M`) or the Moon (`M`) |
| `fx` | `f` | third body true anomaly |
| `fx2` | `F` | partial expression of the third body long-period periodic contribution |
| `fx3` | `F` | partial expression of the third body long-period periodic contribution |
| `third_body_delta_eccentricity` | `e` | long-period periodic contribution of the Sun (`e`) or the Moon (`e`) to the eccentricity |
| `third_body_delta_inclination` | `I` | long-period periodic contribution of the Sun (`I`) or the Moon (`I`) to the inclination |
| `third_body_delta_mean_mootion` | `M` | long-period periodic contribution of the Sun (`M`) or the Moon (`M`) to the mean motion |
| `px4` | `p` | partial expression of the long-period periodic contribution of the Sun (`p`) or the Moon (`p`) to the right ascension of the ascending node and the argument of perigee |
| `px5` | `p` | partial expression of the long-period periodic contribution of the Sun (`p`) or the Moon (`p`) to the right ascension of the ascending node |
### Mathematical expressions
1. [UT1 to Julian conversion](#ut1-to-julian-conversion)
2. [Common initialization](#common-initialization)
3. [Near earth initialization](#near-earth-initialization)
4. [High altitude near earth initialization](#high-altitude-near-earth-initialization)
5. [Elliptic high altitude near earth initialization](#elliptic-high-altitude-near-earth-initialization)
6. [Deep space initialization](#deep-space-initialization)
7. [Third body perturbations](#third-body-perturbations)
8. [Resonant deep space initialization](#resonant-deep-space-initialization)
9. [Geosynchronous deep space initialization](#geosynchronous-deep-space-initialization)
10. [Molniya deep space initialization](#molniya-deep-space-initialization)
11. [Common propagation](#common-propagation)
12. [Near earth propagation](#near-earth-propagation)
13. [High altitude near earth propagation](#high-altitude-near-earth-propagation)
14. [Deep space propagation](#deep-space-propagation)
15. [Third body propagation](#third-body-propagation)
16. [Resonant deep space propagation](#resonant-deep-space-propagation)
17. [Lyddane deep space propagation](#lyddane-deep-space-propagation)
---
#### UT1 to Julian conversion
The epoch (Julian years since UTC 1 January 2000 12h00) can be calculated with either the AFSPC formula:
```
y = (367 y - 7 (y + (m + 9) / 12) / 4 + 275 m / 9 + d
+ 1721013.5
+ (((ns / 10 + s) / 60 + min) / 60 + h) / 24
- 2451545)
/ 365.25
```
or the more accurate version of the same formula:
```
y = (367 y - 7 (y + (m + 9) / 12) / 4 + 275 m / 9 + d - 730531) / 365.25
+ (3600 h + 60 min + s - 43200) / (24 60 60 365.25)
+ ns / (24 60 60 365.25 10)
```
#### Common initialization
```
a = (k / n)
3 3 cosI
p = - J -----------
4 (1 e)
= p / a
= p / (a (1 - / - - / ))
n" = n / (1 + )
p = cos I
p = 1 e
k = 3 p - 1
a" = (k / n")
p = a" (1 - e)
p = a (p - 1)
p = 20 if p < 98
p - 78 if 98 p < 156
78 otherwise
s = p / a + 1
p = ((120 - p) / a)
= 1 / (a" - s)
p = p
= a" e
p = |1 - |
p = p / p
C = B* p n" (a" (1 + / + e (4 + ))
+ / J k (8 + 3 (8 + )) / p)
p = (a" p)
= p
p = / J p n"
p = / p J p
p = - / J p n"
p = - p p + (/ p (4 - 19 p) + 2 p (3 - 7 p)) p
k = - / p (1 - 5 p) + / p (7 - 114 p + 395 p)
p = n" + / p k + / p (13 - 78 p + 137 p)
C = 2 B* n" p a" p (
(2 + / )
+ e (/ + 2 )
- J / (a p) (-3 k (1 - 2 e + (/ - / e ))
+ / (1 - p) (2 - e (1 + )) cos 2 )
k = - / p p p C
k = / C
= p if n" > 2 / 225
p + ( + ) otherwise
= k if n" > 2 / 225
k + ( + ) otherwise
= p if n" > 2 / 225
p + ( + ) otherwise
```
#### Near earth initialization
Defined only if `n" > 2 / 225` (near earth).
```
1 J
k = - - -- sin I
2 J
k = 1 - p
k = 7 p - 1
1 J 3 + 5 p
k = - - -- sin I -------- if |1 + p| > 1.5 10
4 J 1 + p
1 J 3 + 5 p
- - -- sin I ----------- otherwise
4 J 1.5 10
```
#### High altitude near earth initialization
Defined only if `n" > 2 / 225` (near earth) and `p 220 / (a + 1)` (high altitude).
```
D = 4 a" C
p = D C / 3
D = (17 a + s) p
D = / p a" (221 a" + 31 s) C
C = 2 B* p a" p (1 + 2.75 ( + e) + e )
k = (1 + cos M)
k = sin M
k = D + 2 C
k = / (3 D + C (12 D + 10 C))
k = / (3 D + 12 C D + 6 D + 15 C (2 D + C))
```
#### Elliptic high altitude near earth initialization
Defined only if `n" > 2 / 225` (near earth), `p 220 / (a + 1)` (high altitude) and `e > 10` (elliptic).
```
J p n" sin I
k = - 2 B* cos -- ----------------
J e
2 p B*
k = - - -----
3 e
```
#### Deep space initialization
Defined only if `n" 2 / 225` (deep space).
```
e = 365.25 (t + 100)
sin I = 0.39785416
cos I = 0.91744867
sin( - ) = sin
cos( - ) = cos
sin = -0.98088458
cos = 0.1945905
M = (6.2565837 + 0.017201977 e) rem 2
= 4.523602 - 9.2422029 10 e rem 2
cos I = 0.91375164 - 0.03568096
sin I = (1 - cosI)
sin = 0.089683511 sin / sin I
cos = (1 - sin)
= 5.8351514 + 0.001944368 e
0.39785416 sin / sin I
+ tan ------------------------------------------ -
cos cos + 0.91744867 sin sin
sin( - ) = sin cos - cos sin
cos( - ) = cos cos + sin sin
M = (-1.1151842 + 0.228027132 e) rem 2
```
#### Third body perturbations
Defined only if `n" 2 / 225` (deep space).
The following variables are evaluated for two third bodies, the Sun (solar perturbations `s`) and the Moon (lunar perturbations `l`). Variables specific to the third body are annotated with `x`. In other sections, `x` is either `s` or `l`.
```
a = cos cos( - ) + sin cos I sin( - )
a = - sin cos( - ) + cos cos I sin( - )
a = - cos sin( - ) + sin cos I cos( - )
a = sin sin I
a = sin sin( - ) + cos cos I cos( - )
a = cos sin I
a = a cos i + a sin I
a = a cos i + a sin I
a = - a sin I + a cos I
a = - a sin I + a cos I
X = a cos + a sin
X = a cos + a sin
X = - a sin + a cos
X = - a sin + a cos
X = a sin
X = a sin
X = a cos
X = a cos
Z = 12 X - 3 X
Z = 24 X X - 6 X X
Z = 12 X - 3 X
Z = - 6 a a + e (- 24 X X - 6 X X)
Z = - 6 a a + e (-24 X X - 6 X X)
Z = 6 a a + e (24.0 X X - 6 X X)
Z = 6 a a + e (24 X X - 6 X X)
Z = 2 (3 (a + a) + Z e) + p Z
Z = 2 (3 (a + a) + Z e) + p Z
p = C / n"
1 p
p = - - ---
2
p = p
p = - 15 e p
= 0 if I < 5.2359877 10
or I > - 5.2359877 10
- n p (Z + Z) / sin I otherwise
k = 2 p (X X + X X)
k = 2 p (X X - X X)
k = 2 p (- 6 (a a + a a) + e (- 24 (X X + X X) - 6 (X X + X X)))
k = 2 p (Z - Z)
k = - 2 p (2 (6 (a a + a a) + Z e) + p Z)
k = - 2 p (Z - Z)
k = - 2 p (- 21 - 9 e) e
k = 2 p Z
k = 2 p (Z - Z)
k = - 18 p e
k = - 2 p (6 (a a + a a) + e (24 (X X + X X) - 6 (X X + X X)))
k = - 2 p (Z - Z)
I = p n (Z + Z)
= p n (X X + X X)
= p n (Z + Z - 6) - cos I
= - n p (Z + Z - 14 - 6 e)
```
#### Resonant deep space initialization
Defined only if `n" 2 / 225` (deep space) and either:
- `0.0034906585 < n" < 0.0052359877` (geosynchronous)
- `8.26 10 n" 9.24 10` and `e 0.5` (Molniya)
The sidereal time `` at epoch can be calculated with either the IAU formula:
```
c = y / 100
= / ( / 180) (- 6.2 10 c + 0.093104 c
+ (876600 3600 + 8640184.812866) c + 67310.54841) mod 2
```
or the AFSPC formula:
```
d = 365.25 (y + 30) + 1
= 1.7321343856509374 + 1.72027916940703639 10 d + 10
+ (1.72027916940703639 10 + 2) (d - d + 10)
+ 5.07551419432269442 10 d mod 2
```
```
= M + + rem 2 if geosynchronous
M + 2 2 rem 2 otherwise
= p + (k + p) + ( + ) + ( + ) + ( + ) - n" if geosynchronous
p + ( + ) + 2 (p + ( + ) - ) - n" otherwise
```
#### Geosynchronous deep space initialization
Defined only if `n" 2 / 225` (deep space) and `0.0034906585 < n" < 0.0052359877` (geosynchronous orbit).
```
p = 3 (n / a")
= p (/ sinI (1 + 3 p) - / (1 + p))
(1 + 2 e) 2.1460748 10 / a"
= 2 p (/ (1 + p))
(1 + e (- / + / e)) 1.7891679 10
= 3 p (/ (1 + p)) (1 + e (- 6 + 6.60937 e))
2.2123015 10 / a"
```
#### Molniya deep space initialization
Defined only if `n" 2 / 225` (deep space) and `8.26 10 n" 9.24 10` and `e 0.5` (Molniya).
```
p = 3 n" / a"
p = p / a"
p = p / a"
p = p / a"
F = / (1 + 2 p + p)
G = 3.616 - 13.247 e + 16.29 e if e 0.65
- 72.099 + 331.819 e - 508.738 e + 266.724 e otherwise
G = - 19.302 + 117.39 e - 228.419 e + 156.591 e if e 0.65
- 346.844 + 1582.851 e - 2415.925 e + 1246.113 e otherwise
G = - 18.9068 + 109.7927 e - 214.6334 e + 146.5816 e if e 0.65
- 342.585 + 1554.908 e - 2366.899 e + 1215.972 e otherwise
G = - 41.122 + 242.694 e - 471.094 e + 313.953 e if e 0.65
- 1052.797 + 4758.686 e - 7193.992 e + 3651.957 e otherwise
G = - 146.407 + 841.88 e - 1629.014 e + 1083.435 e if e 0.65
- 3581.69 + 16178.11 e - 24462.77 e + 12422.52 e otherwise
G = - 532.114 + 3017.977 e - 5740.032 e + 3708.276 e if e 0.65
1464.74 - 4664.75 e + 3763.64 e if 0.65 < e < 0.715
- 5149.66 + 29936.92 e - 54087.36 e + 31324.56 e otherwise
G = - 853.666 + 4690.25 e - 8624.77 e + 5341.4 e if e < 0.7
- 40023.88 + 170470.89 e - 242699.48 e + 115605.82 e otherwise
G = - 822.71072 + 4568.6173 e - 8491.4146 e + 5337.524 e if e < 0.7
- 51752.104 + 218913.95 e - 309468.16 e + 146349.42 e otherwise
G = - 919.2277 + 4988.61 e - 9064.77 e + 5542.21 e if e < 0.7
- 37995.78 + 161616.52 e - 229838.2 e + 109377.94 e otherwise
D = p 1.7891679 10 F (- 0.306 - 0.44 (e - 0.64))
D = p 1.7891679 10 (/ sinI) G
D = p 3.7393792 10 (/ sin I (1 - 2 p - 3 p)) G
D = p 3.7393792 10 (- / sin I (1 + 2 p - 3 p)) G
D = 2 p 7.3636953 10 (35 sinI F) G
D = 2 p 7.3636953 10 (/ sinI) G
D = p 1.1428639 10 (/ sin I
(sinI (1 - 2 p - 5 p)
+ 0.33333333 (- 2 + 4 p + 6 p))) G
D = p 1.1428639 10 (sin I
(4.92187512 sinI (- 2 - 4 p + 10 p)
+ 6.56250012 (1 + p - 3 p))) G
D = 2 p 2.1765803 10 (/ sin I
(2 - 8 p + p (- 12 + 8 p + 10 p))) G
D = 2 p 2.1765803 10 (/ sin I
(- 2 - 8 p + p (12 + 8 p - 10 p))) G
```
#### Common propagation
The following values depend on the propagation time `t` (minutes since epoch).
Named conditions have the following meaning:
- `near earth`: `n" 2 / 225`
- `low altitude near earth`: `near earth` and `p < 220 / (a + 1)`
- `high altitude near earth`: `near earth` and `p 220 / (a + 1)`
- `elliptic high altitude near earth`: `high altitude near earth` and `e > 10`
- `non-elliptic near earth`: `low altitude near earth` or `high altitude near earth` and `e 10`
- `deep space`: `n" > 2 / 225`
- `non-Lyddane deep space`: `deep space` and `I 0.2`
- `Lyddane deep space`: `deep space` and `I < 0.2`
- `AFSPC Lyddane deep space`: `Lyddane deep space` and use the same expression as the original AFSPC implementation, with an `` discontinuity at `p = 0`
```
p = + t + k t
p = + t
I = I if near earth
I + I t + (I + I) otherwise
= p if near earth
p + (p + p) / sin I if non-Lyddane deep space
p + 2 if Lyddane deep space and p + < p rem 2
p - 2 if Lyddane deep space and p - > p rem 2
p otherwise
e = 10 if near earth and p < 10
p if near earth and p 10
10 + (e + e) if deep space and p < 10
p + (e + e) otherwise
= p - p if elliptic high altitude near earth
p if non-elliptic near earth
p + (p + p) - cos I (p + p) / sin I if non-Lyddane deep space
p + (p + p) + cos I ((p rem 2) - )
- (I + I) (p mod 2) sin I if AFSPC Lyddane deep space
p + (p + p) + cos I ((p rem 2) - )
- (I + I) (p rem 2) sin I otherwise
M = p + n" (k t + k t + t (k + t k) if high altitude near earth
p + n" k t if low altitude near earth
p + (M + M) + n" k t otherwise
a = a" (1 - C t - D t - D t - D t) if high altitude near earth
a" (1 - C t) if low altitude near earth
p (1 - C t) otherwise
n = k / a
p = k if near earth
1 J
- - -- sin I othewise
2 J
p = k if near earth
1 - cosI othewise
p = k if near earth
7 cosI - 1 otherwise
p = k if near earth
1 J 3 + 5 cos I
- - -- sin I ----------- if deep space and |1 + cos I| > 1.5 10
4 J 1 + cos I
1 J 3 + 5 cos I
- - -- sin I ----------- otherwise
4 J 1.5 10
p = k if near earth
3 cosI - 1 otherwise
p = 1 / (a (1 - e))
a = e cos
a = e sin + p p
p = M + + p p a rem 2
(E + ) = p
p - a cos (E + ) + a sin (E + ) - (E + )
(E + ) = ---------------------------------------------------
1 - cos (E + ) a - sin (E + ) a
(E + ) = (E + ) + (E + )|[-0.95, 0.95]
E + = (E + ) if j [0, 9], (E + ) 10
(E + ) otherwise, with j the smallest integer | (E + ) < 10
p = a + a
p = a (1 - p)
p = a sin(E + ) - a cos(E + )
r = a (1 - (a cos(E + ) + a sin(E + )))
r = a p / r
= (1 - p)
p = p / (1 + )
p = a / r (sin(E + ) - a - a p)
p = a / r (cos(E + ) - a + a p)
p
u = tan ---
p
p = 2 p p
p = 1 - 2 p
p = (/ J / p) / p
r = r (1 - / p p) + / (/ J / p) p p
u = u - / p p p
= + / p cos I p
I = I + / p cos I sin I p
r = r + n (/ J / p) p / k
rf = p / r + n (/ J / p) (p p + / p) / k
u = - sin cos I sin u + cos cos u
u = cos cos I sin u + sin cos u
u = sin I sin u
r = r u a
r = r u a
r = r u a
r = (r u + rf (- sin cos I cos u - cos sin u)) a k / 60
r = (r u + rf (cos cos I cos u - sin sin u)) a k / 60
r = (r u + rf (sin I cos u)) a k / 60
```
#### Near earth propagation
Defined only if `n" > 2 / 225` (near earth).
```
p = M + t
p = | e - (C t + C (sin p - k)) if high altitude
| e - C t otherwise
```
#### High altitude near earth propagation
Defined only if `n" > 2 / 225` (near earth) and `p 220 / (a + 1)` (high altitude).
`elliptic` means `e > 10`.
```
p = k ((1 + cos p) - k) + k t
p = p + p if elliptic
p otherwise
```
#### Deep space propagation
Defined only if `n" 2 / 225` (deep space).
```
p = (k / (n + n (t - t) + / n (t - t))) if geosynchronous or Molniya
a" otherwise
p = + (t - t) + / n (t - t) - p - p + if geosynchronous
+ (t - t) + / n (t - t) - 2 p + 2 if Molniya
M + t otherwise
j is the largest positive integer | t t if t > 0
the smallest negative integer | t t if t < 0
0 otherwise
p = e + t - C t
```
#### Third body propagation
Defined only if `n" 2 / 225` (deep space).
The following variables are evaluated for two third bodies, the Sun (solar perturbations `s`) and the Moon (lunar perturbations `l`). Variables specific to the third body are annotated with `x`. In other sections, `x` is either `s` or `l`.
```
M = M + n t
f = M + 2 e sin M
F = / sinf - /
F = - / sin f cos f
e = k F + k F
I = k F + k F
M = k F + k F + k sin f
p = k F + k F + k sin f
p = k F + k F
```
#### Resonant deep space propagation
Defined only if `n" 2 / 225` (deep space) and either:
- `0.0034906585 < n" < 0.0052359877` (geosynchronous)
- `8.26 10 n" 9.24 10` and `e 0.5` (Molniya)
```
= + 4.37526908801129966 10 t rem 2
t = |t| if t > 0
-|t| if t < 0
0 otherwise
= n +
n = sin( - ) + sin(2 ( - )) + sin(3 ( - )) if geosynchronous
D sin((l - 2 p) + m / 2 - G) otherwise
n = ( cos( - ) + cos(2 ( - )) + cos(3 ( - ))) if geosynchronous
( m / 2 D cos((l - 2 p) + m / 2 - G)) otherwise
(l, m, p, k) {(2, 2, 0, -1), (2, 2, 1, 1), (3, 2, 1, 0),
(3, 2, 2, 2), (4, 4, 1, 0), (4, 4, 2, 2), (5, 2, 2, 0),
(5, 2, 3, 2), (5, 4, 2, 1), (5, 4, 3, 3)}
t = t + t
n = n + n t + n (t / 2)
= + t + n (t / 2)
```
#### Lyddane deep space propagation
Defined only if `n" 2 / 225` (deep space) and `I < 0.2` (Lyddane).
```
sin I sin p + (p + p) cos p + (I + I) cos I sin p
p = tan -------------------------------------------------------------
sin I cos p - (p + p) sin p + (I + I) cos I cos p
```
## References
[1] David A. Vallado, Paul Crawford, R. S. Hujsak and T. S. Kelso, "Revisiting Spacetrack Report #3", presented at the AIAA/AAS Astrodynamics Specialist Conference, Keystone, CO, 2006 August 2124, https://doi.org/10.2514/6.2006-6753
[2] F. R. Hoots, P. W. Schumacher Jr. and R. A. Glover, "History of Analytical Orbit Modeling in the U. S. Space Surveillance System", Journal of Guidance, Control, and Dynamics, 27(2), 174185, 2004, https://doi.org/10.2514/1.9161
[3] F. R. Hoots and R. L. Roehrich, "Spacetrack Report No. 3: Models for propagation of NORAD element sets", U.S. Air Force Aerospace Defense Command, Colorado Springs, CO, 1980, https://www.celestrak.com/NORAD/documentation/
[4] R. S. Hujsak, "A Restricted Four Body Solution for Resonating Satellites Without Drag", Project SPACETRACK, Rept. 1, U.S. Air Force Aerospace Defense Command, Colorado Springs, CO, Nov. 1979, https://doi.org/10.21236/ada081263
Owner
- Name: Virgile Andreani
- Login: Armavica
- Kind: user
- Location: Boston
- Repositories: 60
- Profile: https://github.com/Armavica
Hi! I am a physicist by training with a PhD in computational biology. I love Rust, Monte-Carlo methods and scientific computing.