SGP4 — Rust 实现版本 // Lib.rs • Rust 包仓库

22 个版本 (6 个稳定版本)

2.2.0	2024 年 5 月 25 日
2.1.0	2023 年 11 月 16 日
1.2.2	2023 年 3 月 5 日
0.9.1	2023 年 1 月 30 日
0.3.1	2020 年 7 月 16 日

#40 in 算法

715 次月下载
用于 3 个crate（2 个直接使用）

MIT 许可证

250KB
2.5K SLoC

SGP4 算法，从 Celestrak 的参考实现移植到 Rust。

代码完全重构，以利用 Rust 的代数数据类型，并突出实现与参考数学方程之间的关系。

SGP4 可以通过 WebAssembly 包装器从 JavaScript 或 Python 调用。有关如何安装和使用 SGP4 作为 WAPM 包的信息，请参阅 https://github.com/wasmerio/sgp4。

数值预测几乎与 Celestrak 的相同。观察到的差异（在历元后三年半，位置小于 2 × 10⁻⁷ km，速度小于 10⁻⁹ km.s⁻¹）远低于算法的精度。

我们从不完整的 https://github.com/natronics/rust-sgp4 吸取灵感，使用 UTF-8 字符编写数学表达式。

文档
无 std 或 alloc 的环境
基准测试
变量和数学表达式
- 变量
- 数学表达式
参考文献

文档

代码文档托管在 https://docs.rs/sgp4/latest/sgp4/。

示例可以在本仓库的 examples 目录中找到

examples/celestrak.rs 从 Celestrak 获取最新的伽利略 OMMs 并进行传播
examples/omm.rs 解析并传播 JSON 编码的 OMM
examples/space-track.rs 从 Space-Track 获取最近的 20 次发射 OMMs 并进行传播
examples/tle.rs 解析并传播 TLE
examples/tle_afspc.rs 使用 AFSPC 兼容模式解析和传播 TLE
examples/advanced.rs 利用高级 API 以（略微）加速深空共振卫星的传播

要运行示例（此处为 examples/celestrak.rs），使用

cargo run --example celestrak

要运行 Space-Track 示例，您必须首先将 Space-Track.org 凭证分配给字段 identity 和 password（参见 examples/space-track.rs 中的第 3 和 4 行）。

无 std 或 alloc 的环境

此 crate 支持 no_std 环境。TLE 解析和 SGP4 传播也不需要 alloc。当 std 不可用时，我们使用 num-traits 与 libm 一起使用浮点函数。

所有与 serde 相关的功能，尤其是 OMM 解析，都需要 alloc。

基准测试

基准代码可在 https://github.com/neuromorphicsystems/sgp4-benchmark 找到。它比较了不同配置下的两个 SGP4 实现

cpp：改进模式下的 Celestrak 实现 [1]
cpp-afspc：AFSPC 兼容模式下的 Celestrak 实现 [1]
cpp-fastmath：带有 fast-math 编译器标志的改进模式下的 Celestrak 实现 [1]
cpp-afspc-fastmath：带有 fast-math 编译器标志的 AFSPC 兼容模式下的 Celestrak 实现 [1]
rust：默认模式下的我们的 Rust 实现
rust-afspc：AFSPC 兼容模式下的我们的 Rust 实现

此基准测试不应与本仓库 bench 目录中的代码混淆。后者仅考虑了 Celestrak 目录的一小部分（[1] 中推荐测试）并且没有测量原始的 C++ 实现。

以下为本机配置

CPU - Intel Core i7-8700 @ 3.20GHz
RAM - Kingston DDR4 @ 2.667 GHz
OS - Ubuntu 16.04
编译器 - Rust 1.44.1 和 gcc 9.3.0

精度衡量的是每个实现相对于参考实现（cpp-afspc）在完整的 Celestrak 目录中的最大传播误差（24小时内的 1 分钟时间步长）。

实现	最大位置误差	最大速度误差
`cpp-afspc`	（参考）	（参考）
`cpp`	1.05 公里	1.30 × 10⁻³ 公里·秒⁻¹
`cpp-fastmath`	1.05 公里	1.30 × 10⁻³ 公里·秒⁻¹
`cpp-afspc-fastmath`	4.21 × 10⁻⁸ 公里	7.51 × 10⁻¹² 公里·秒⁻¹
`rust`	1.05 公里	1.30 × 10⁻³ 公里·秒⁻¹
`rust-afspc`	4.19 × 10⁻⁸ 公里	7.46 × 10⁻¹² 公里·秒⁻¹

Rust 和 C++ 快速数学误差具有相同数量级。在两种情况下，它们都可以归因于使用不同的浮点运算实现数学上相同的表达式。

速度衡量的是使用单个线程传播 Celestrak 目录中的每个卫星所需的时间（24小时内的 1 分钟时间步长）。每个实现采样 100 个值。

实现	最小	Q1	中位数	Q3	最大	相对差异
`cpp-afspc`	8.95 秒	9.02 秒	9.03秒	9.06秒	9.18秒	（参考）
`cpp`	8.95 秒	9.01秒	9.04秒	9.06秒	9.25秒	+ 0 %
`cpp-fastmath`	7.67秒	7.74秒	7.77秒	7.79秒	7.90秒	- 14 %
`cpp-afspc-fastmath`	7.70秒	7.74秒	7.76秒	7.79秒	7.86秒	- 14 %
`rust`	8.36秒	8.41秒	8.43秒	8.45秒	8.53秒	- 7 %
`rust-afspc`	8.36秒	8.41秒	8.43秒	8.46秒	8.59秒	- 7 %

Rust快速数学支持仍在开发中（参见https://github.com/rust-lang/rust/issues/21690）。与C++类似，它应该对精度的影响非常小，同时提供显著的速度提升。

变量和数学表达式

变量

每个变量都用于存储一个表达式的结果，且仅存储一个表达式的结果。大多数变量是不可变的，例外的是用于解开普勒方程的变量(E + ω)ᵢ以及用于积分地球重力共振效应的状态变量tᵢ、nᵢ和λᵢ。

以下表格列出了代码中使用的变量及其相关的数学符号。尽可能使用来自[2]的符号。在[2]中未命名的部分表达式，如果是初始化和传播中共享的，则遵循惯例kₙ, n ∈ ℕ；如果是初始化或传播局部的，则遵循惯例pₙ, n ∈ ℕ。

初始化变量
传播变量
第三体初始化变量
第三体传播变量

初始化变量

以下变量完全取决于历元元素。

变量	符号	描述
`元素::时间.年()`	`yᵤ`	格里高利历年份
`元素::时间.月()`	`mᵤ`	格里高利历月份，范围在`[1, 12]`
`元素::时间.日()`	`dᵤ`	格里高利历日，范围在`[1, 31]`
`元素::时间.时()`	`hᵤ`	午夜过后的小时数，范围在`[0, 23]`
`元素::时间.分()`	`minᵤ`	自小时起经过的分钟数，范围在`[0, 59]`
`元素::时间.秒()`	`sᵤ`	自分钟起经过的秒数，范围在`[0, 59]`
`元素::时间.纳秒()`	`nsᵤ`	自秒起经过的纳秒数，范围在`[0, 10⁹[`
`历元`	`y₂₀₀₀`	自UTC 1月1日2000年12点00分以来的儒略年（J2000）
`d1900`	`d₁₉₀₀`	自UTC 1月1日1900年12点00分以来的儒略日（J1900）
`d1970`	`d₁₉₇₀`	自UTC 1月1日1970年12点00分以来的儒略日（J1970）
`c2000`	`c₂₀₀₀`	自UTC 1月1日2000年12点00分以来的儒略世纪（J2000）
`位能.ae`	`aₑ`	地球赤道半径，单位为千米
`位能.ke`	`kₑ`	地球引力参数的平方根，单位为地球半径³ min⁻²
`位能.j2`	`J₂`	未归一化的第二带谐系数
`位能.j3`	`J₃`	未归一化的第三带谐系数
`位能.j4`	`J₄`	未归一化的第四区谐波
`kozai_mean_motion`	`n₀`	每日常数（Kozai convention）在历元时的弧度/分钟⁻¹
`a1`	`a₁`	历元的半长轴（Kozai convention）
`p0`	`p₀`	𝛿₀ 和 𝛿₁ 的部分表达式
`d1`	`𝛿₁`	用于 Kozai 到 Brouwer 转换
`d0`	`𝛿₀`	用于 Kozai 到 Brouwer 转换
`B*`	`B*`	辐射压力系数（地球半径⁻¹）
`orbit_0.倾角`	`I₀`	赤道与轨道平面在历元时的夹角（弧度）
`orbit_0.right_ascension`	`Ω₀`	春分点与轨道穿过赤道平面的点在历元时的夹角（弧度）
`orbit_0.eccentricity`	`e₀`	历元的轨道形状
`orbit_0.argument_of_perigee`	`ω₀`	升交点与地球最近点的夹角在历元时的弧度
`orbit_0.mean_anomaly`	`M₀`	卫星位置从近地点测量的角度在历元时（弧度）
`orbit_0.mean_motion`	`n₀"`	每日常数（Brouwer convention）在历元时的弧度/分钟⁻¹
`p1`	`p₁`	历元倾角的余弦值，在初始化过程中用于多个表达式（[2]中的θ，为了避免与恒星时混淆而更名）
`p2`	`p₂`	多个初始化表达式的部分表达式
`a0`	`a₀"`	历元的半长轴（Brouwer convention）
`p3`	`p₃`	近地点（地球半径）
`p4`	`p₄`	近地点高度（千米）
`p5`	`p₅`	关于 `s` 的部分表达式
`s`	`s`	大气阻力的轨道高度参数
`p6`	`p₆`	大气阻力的部分表达式
`xi`	`ξ`	多个初始化表达式的部分表达式
`p7`	`p₇`	多个初始化表达式的部分表达式
`eta`	`η`	多个初始化表达式、近地点引数和偏心率在近地高空的传播中的部分表达式
`p8`	`p₈`	多个初始化表达式的部分表达式
`p9`	`p₉`	多个初始化表达式的部分表达式
`c1`	`C₁`	多个初始化和传播表达式的部分表达式
`p10`	`p₁₀`	多个初始化表达式的部分表达式
`b0`	`β₀`	多个初始化表达式的部分表达式
`p11`	`p₁₁`	多个初始化表达式的部分表达式
`p12`	`p₁₂`	多个初始化表达式的部分表达式
`p13`	`p₁₃`	多个初始化表达式的部分表达式
`p14`	`p₁₄`	多个初始化表达式的部分表达式
`p15`	`p₁₅`	多个初始化表达式的部分表达式
`k14`	`k₁₄`	添加太阳和月球摄动前的近地点引数的一阶系数
`c4`	`C₄`	多个初始化和传播表达式的部分表达式（与[2]中的C₄表达式不同，乘以因子B*）
`right_ascension_dot`	`Ω̇`	右升交点的一阶系数
`argument_of_perigee_dot`	`ω̇`	近地点引数的一阶系数
`mean_anomaly_dot`	`Ṁ`	偏心率的一阶系数
`k0`	`k₀`	添加摄动前的右升交点的二阶系数
`k1`	`k₁`	偏心率二阶系数的部分表达式
`k2`	`k₂`	近地传播中关于 `aᵧₙ` 的部分表达式
`k3`	`k₃`	近地传播中关于 `rₖ`、`ṙₖ` 和 `rḟₖ` 的部分表达式
`k4`	`k₄`	近地传播中关于 `uₖ` 的部分表达式
`k5`	`k₅`	近地传播中关于初始开普勒变量 `p₃₈` 的部分表达式
`k6`	`k₆`	多个初始化表达式和在近地传播中关于 `rₖ` 和 `rḟₖ` 的部分表达式
`第2个变量d`	`D₂`	多个近地轨道初始化表达式和半长轴在近地传播中的部分表达式
`p16`	`p₁₆`	多个近地轨道初始化表达式的部分表达式
`d3`	`D₃`	多个近地轨道初始化表达式和半长轴在近地传播中的部分表达式
`d4`	`D₄`	多个近地轨道初始化表达式和半长轴在近地传播中的部分表达式
`c5`	`C₅`	多个初始化和传播表达式的部分表达式（与[2]中的`C₅`表达式相差一个因子B*）
`k7`	`k₇`	历元平均经度的正弦
`k8`	`k₈`	高空近地轨道传播中平均经度三阶系数的部分表达式
`k9`	`k₉`	高空近地轨道传播中平均经度四阶系数的部分表达式
`k10`	`k₁₀`	高空近地轨道传播中平均经度五阶系数的部分表达式
`k11`	`k₁₁`	在离心高空近地轨道传播中近地点和平均经度的部分表达式
`k12`	`k₁₂`	在离心高空近地轨道传播中近地点和平均经度的部分表达式
`k13`	`k₁₃`	在离心高空近地轨道传播中近地点和平均经度的部分表达式
`月球赤经ε`	`Ωₗₑ`	月球升交点赤经
`月球赤经正弦`	`sin Ωₗ`	以赤道为参考的月球升交点赤经的正弦
`月球赤经余弦`	`cos Ωₗ`	以赤道为参考的月球升交点赤经的余弦
`月球近地点经度`	`ωₗ`	月球近地点经度
`历元0的恒星时`	`θ₀`	历元的格林尼治恒星时
`lambda_0`	`λ₀`	历元的地球重力共振变量
`lambda_dot_0`	`λ̇₀`	历元地球重力共振变量的时间导数
`p17`	`p₁₇`	`𝛿ᵣ₁`、`𝛿ᵣ₂`和`𝛿ᵣ₃`的部分表达式
`dr1`	`𝛿ᵣ₁`	地球同步轨道的第一个地球重力共振系数（[2]中的`𝛿₁`，重命名为避免与Kozai到Brouwer转换中使用的`𝛿₁`混淆）
`dr2`	`𝛿ᵣ₂`	地球同步轨道的第二个地球重力共振系数（[2]中的`𝛿₂`，重命名为匹配`𝛿ᵣ₁`）
`dr3`	`𝛿ᵣ₃`	地球同步轨道的第三个地球重力共振系数（[2]中的`𝛿₃`，重命名为匹配`𝛿ᵣ₁`）
`p18`	`p₁₈`	`D₂₂₀₋₁`和`D₂₂₁₁`的部分表达式
`p19`	`p₁₉`	`D₃₂₁₀`和`D₃₂₂₂`的部分表达式
`p20`	`p₂₀`	`D₄₄₁₀`和`D₄₄₂₂`的部分表达式
`p21`	`p₂₁`	`D₅₂₂₀`、`D₅₂₃₂`、`D₅₄₂₁`和`D₅₄₃₃`的部分表达式
`f220`	`F₂₂₀`	`D₂₂₀₋₁`和`D₄₄₁₀`的部分表达式
`g211`	`G₂₁₁`	`D₂₂₁₁`的部分表达式
`g310`	`G₃₁₀`	`D₃₂₁₀`的部分表达式
`g322`	`G₃₂₂`	`D₃₂₂₂`的部分表达式
`g410`	`G₄₁₀`	`D₄₄₁₀`的部分表达式
`g422`	`G₄₂₂`	`D₄₄₂₂`的部分表达式
`g520`	`G₅₂₀`	`D₅₂₂₀`的部分表达式
`g532`	`G₅₃₂`	表达式部分 `D₅₂₃₂`
`g521`	`G₅₂₁`	表达式部分 `D₅₄₂₁`
`g533`	`G₅₃₃`	表达式部分 `D₅₄₃₃`
`d220-1`	`D₂₂₀₋₁`	对于莫尼娅轨道的重力共振系数（[2]中的`Dₗₘₚₖ`表达式缺少来自[4]原始方程的系数 `l - 2p + k`，其中 `k = -1` 而不是 `1`）
`d2211`	`D₂₂₁₁`	对于莫尼娅轨道的重力共振系数（[2]中的`Dₗₘₚₖ`表达式缺少来自[4]原始方程的系数 `l - 2p + k`）
`d3210`	`D₃₂₁₀`	参见 `D₂₂₁₁`
`d3222`	`D₃₂₂₂`	参见 `D₂₂₁₁`
`d4410`	`D₄₄₁₀`	参见 `D₂₂₁₁`
`d4422`	`D₄₄₂₂`	参见 `D₂₂₁₁`
`d5220`	`D₅₂₂₀`	参见 `D₂₂₁₁`
`d5232`	`D₅₂₃₂`	参见 `D₂₂₁₁`
`d5421`	`D₅₄₂₁`	参见 `D₂₂₁₁`
`d5433`	`D₅₄₃₃`	参见 `D₂₂₁₁`

传播变量

以下表达式依赖于传播时间 t。

变量	符号	描述
`t`	`t`	自历元起经过的分钟（可以是负数）
`p22`	`p₂₂`	没有地球重力共振、太阳和月亮贡献的升交点赤经
`p23`	`p₂₃`	没有高空气阻效应、地球重力共振和太阳月亮贡献的近地点幅角
`轨道.倾角`	`I`	历元时的倾角加上 `t`，没有地球重力短期效应的影响
`轨道.right_ascension`	`Ω`	历元时的升交点赤经加上 `t`，没有地球重力短期效应的影响
`轨道.eccentricity`	`e`	历元时的偏心率加上 `t`，没有地球重力短期效应的影响
`轨道.argument_of_perigee`	`ω`	历元时的近地点幅角加上 `t`，没有地球重力短期效应的影响
`轨道.mean_anomaly`	`M`	历元时的平均经度加上 `t`，没有地球重力短期效应的影响
`轨道.mean_motion`	`n`	历元时的平均运动加上 `t`，没有地球重力短期效应的影响
`a`	`a`	半长轴
`p32`	`p₃₂`	`aᵧₙ` 的部分表达式
`p33`	`p₃₃`	`rₖ`、`ṙₖ` 和 `rḟₖ` 的部分表达式
`p34`	`p₃₄`	`uₖ` 的部分表达式
`p35`	`p₃₅`	初始开普勒变量 `p₃₈` 的部分表达式
`p36`	`p₃₆`	`rₖ` 和 `rḟₖ` 的部分表达式
`p37`	`p₃₇`	`aᵧₙ` 和初始开普勒变量 `p₃₈` 的部分表达式
`axn`	`aₓₙ`	投影到节线上的归一化线性偏心率
`ayn`	`aᵧₙ`	投影到节线法线上的归一化线性偏心率
`p38`	`p₃₈`	初始开普勒变量（[2]中的 `U`，为了避免与真近点角加上近地点幅角 `u` 混淆而更名）
`ew`	`(E+ω)ᵢ`	在迭代过程中使用开普勒变量来估计偏心异常 `E`
`增量`	`Δ(E+ω)ᵢ`	迭代 `i` 时对开普勒变量的修正
`第39页`	`第39页`	历元处的偏心率加上 `t`
`pl`	`pₗ`	半通径
`第40页`	`第40页`	投影到半短轴上的归一化线性偏心率
`r`	`r`	半径（到焦点的距离）不考虑地球引力的短周期效应
`r_dot`	`ṙ`	半径的时间导数不考虑地球引力的短周期效应
`b`	`β`	半短轴与半长轴之比
`第41页`	`第41页`	`p₄₂` 和 `p₄₃` 的部分表达式
`p42`	`p₄₂`	`u` 的正弦
`p43`	`p₄₃`	`u` 的余弦
`u`	`u`	真实近点角加上近点赤经不考虑地球引力的短周期效应
`第44页`	`第44页`	`sin(2 u)`，`uₖ`，`Ωₖ` 和 `ṙₖ` 的部分表达式
`第45页`	`第45页`	`cos(2 u)`，`rₖ`，`Iₖ` 和 `rḟₖ` 的部分表达式
`第46页`	`第46页`	`rₖ`，`uₖ`，`Iₖ` 和 `Ωₖ` 的部分表达式
`rk`	`rₖ`	半径（到焦点的距离）
`uk`	`uₖ`	真实近点角加上近点赤经
`inclination_k`	`Iₖ`	历元处的倾角加上 `t`
`right_ascension_k`	`Ωₖ`	历元处的升交点赤经加上 `t`
`rk_dot`	`ṙₖ`	半径的时间导数
`rfk_dot`	`rḟₖ`	半径乘以真实近点角的时间导数
`u0`	`u₀`	位置单位向量的x分量
`u1`	`u₁`	位置单位向量的y分量
`u2`	`u₂`	位置单位向量的z分量
`预测.位置[0]`	`r₀`	位置向量的x分量（公里）（真实赤道，平赤经（TEME）的历元参考框架）
`预测.位置[1]`	`r₁`	位置向量的y分量（公里）（真实赤道，平赤经（TEME）的历元参考框架）
`预测.位置[2]`	`r₂`	位置向量的z分量（公里）（真实赤道，平赤经（TEME）的历元参考框架）
`预测.速度[0]`	`ṙ₀`	速度向量的x分量（公里·秒⁻¹）（真实赤道，平赤经（TEME）的历元参考框架）
`预测.速度[1]`	`ṙ₁`	速度向量的y分量（公里·秒⁻¹）（真实赤道，平赤经（TEME）的历元参考框架）
`预测.速度[2]`	`ṙ₂`	速度向量的z分量（公里·秒⁻¹）（真实赤道，平赤经（TEME）的历元参考框架）
`p24`	`p₂₄`	近地传播中不考虑阻力贡献的平均近点角
`p25`	`p₂₅`	近地传播中考虑椭圆修正且不考虑阻力贡献的平均近点角的部分表达式
`p26`	`p₂₆`	近地传播中考虑椭圆修正且不考虑阻力贡献的平均近点角
`p27`	`p₂₇`	近地传播中的非固定偏心率
`p28`	`p₂₈`	深空传播中考虑共振修正的半长轴
`p29`	`p₂₉`	深空传播中的共振修正平均偏差
`p31`	`p₃₁`	深空传播中的非夹紧偏心率
`恒星时`	`θ`	历元加 `t` 的恒星时
`delta_t`	`Δt`	在地球重力共振效应积分中使用的步长，以分钟为单位（要么是 `720`，要么是 `-720`）
`lambda_dot`	`λ̇ᵢ`	在历元加 `i Δt` 时的地球重力共振效应变量的时间导数
`ni_dot`	`ṅᵢ`	在历元加 `i Δt` 时的平均运动时间导数
`ni_ddot`	`n̈ᵢ`	在历元加 `i Δt` 时的平均运动二阶导数
`ResonanceState::t`	`tᵢ`	地球重力积分器的共振效应时间（`i Δt`）
`ResonanceState::mean_motion`	`nᵢ`	在历元加 `Δt i` 时的平均运动时间导数
`ResonanceState::lambda`	`λᵢ`	在历元加 `i Δt` 时的地球重力共振效应变量
`p30`	`p₃₀`	在Lyddane深空传播中的非归一化 `Ω`

第三体初始化变量

计算太阳和月亮对轨道元素的贡献使用了一套独特的表达式。在 src/third_body.rs 中提供了这些表达式的通用实现。特定于第三体（太阳或月亮）的变量用 x 注释。在其他每个文件中，如果这些变量对应于太阳摄动，则用 s 注释，如果对应于月球摄动，则用 l 注释。

aₓₙ、Xₓₙ、Zₓₙ（《n ∈ ℕ》）、Fₓ₂ 和 Fₓ₃ 变量对应于 [2] 中的 aₙ、Xₙ、Zₙ、F₂ 和 F₃ 变量。添加的 x 突出了对摄动第三体的依赖。

以下变量完全取决于历元元素。

变量	符号	描述
`third_body_inclination_sine`	`sin Iₓ`	太阳（`sin Iₛ`）或月亮（`sin Iₗ`）倾角的正弦
`third_body_inclination_cosine`	`cos Iₓ`	太阳（`cos Iₛ`）或月亮（`cos Iₗ`）倾角的余弦
`delta_right_ascension_sine`	`sin(Ω₀-Ωₓ)`	卫星在历元上升节点赤经与太阳（`sin(Ω₀ - Ωₛ)`）或月亮（`sin(Ω₀ - Ωₗ)`）之间差异的正弦
`delta_right_ascension_cosine`	`cos(Ω₀-Ωₓ)`	卫星在历元上升节点赤经与太阳（`cos(Ω₀ - Ωₛ)`）或月亮（`cos(Ω₀ - Ωₗ)`）之间差异的余弦
`third_body_argument_of_perigee_sine`	`sin ωₓ`	太阳（`sin ωₛ`）或月亮（`sin ωₗ`）近地点角正弦
`third_body_argument_of_perigee_cosine`	`cos ωₓ`	太阳（`sin ωₛ`）或月亮（`cos ωₗ`）近地点角余弦
`third_body_mean_anomaly_0`	`Mₓ₀`	太阳（`Mₛ₀`）或月亮（`Mₗ₀`）在历元的平均偏差
`ax1`	`aₓ₁`	多个 `Xₓₙ` 和 `Zₓₙ` 表达式的部分表达式
`ax3`	`aₓ₃`	多个 `Xₓₙ` 和 `Zₓₙ` 表达式的部分表达式
`ax7`	`aₓ₇`	多个 `aₓ₂` 和 `aₓ₅` 的部分表达式
`ax8`	`aₓ₈`	多个 `aₓ₂` 和 `aₓ₅` 的部分表达式
`ax9`	`aₓ₉`	多个 `aₓ₄` 和 `aₓ₆` 的部分表达式
`ax10`	`aₓ₁₀`	多个 `aₓ₄` 和 `aₓ₆` 的部分表达式
`ax2`	`aₓ₂`	多个 `Xₓₙ` 和 `Zₓₙ` 表达式的部分表达式
`ax4`	`aₓ₄`	多个 `Xₓₙ` 和 `Zₓₙ` 表达式的部分表达式
`ax5`	`aₓ₅`	多个 `Xₓₙ` 和 `Zₓₙ` 表达式的部分表达式
`ax6`	`aₓ₆`	多个 `Xₓₙ` 和 `Zₓₙ` 表达式的部分表达式
`xx1`	`Xₓ₁`	多个 `Zₓₙ` 表达式、`kₓ₀`、`kₓ₁` 和 `ėₓ` 的部分表达式
`xx2`	`Xₓ₂`	多个 `Zₓₙ` 表达式、`kₓ₀`、`kₓ₁` 和 `ėₓ` 的部分表达式
`xx3`	`Xₓ₃`	多个 `Zₓₙ` 表达式、`kₓ₀`、`kₓ₁` 和 `ėₓ` 的部分表达式
`xx4`	`Xₓ₄`	多个 `Zₓₙ` 表达式、`kₓ₀`、`kₓ₁` 和 `ėₓ` 的部分表达式
`xx5`	`Xₓ₅`	多个 `Zₓₙ` 表达式的部分表达式
`xx6`	`Xₓ₆`	多个 `Zₓₙ` 表达式的部分表达式
`xx7`	`Xₓ₇`	多个 `Zₓₙ` 表达式的部分表达式
`xx8`	`Xₓ₈`	多个 `Zₓₙ` 表达式的部分表达式
`zx31`	`Zₓ₃₁`	`Zₓ₃`、`kₓ₈` 和 `ω̇ₓ` 的部分表达式
`zx32`	`Zₓ₃₂`	`Zₓ₂`、`kₓ₇` 和 `ω̇ₓ` 的部分表达式
`zx33`	`Zₓ₃₃`	`Zₓ₃`、`kₓ₈` 和 `ω̇ₓ` 的部分表达式
`zx11`	`Zₓ₁₁`	`kₓ₃` 和 `İₓ` 的部分表达式
`zx13`	`Zₓ₁₃`	`kₓ₃` 和 `İₓ` 的部分表达式
`zx21`	`Zₓ₂₁`	`kₓ₁₁` 和 `Ω̇ₓ` 的部分表达式
`zx23`	`Zₓ₂₃`	`kₓ₁₁` 和 `Ω̇ₓ` 的部分表达式
`zx1`	`Zₓ₁`	`kₓ₅` 和 `Ṁₓ` 的部分表达式
`zx3`	`Zₓ₃`	`kₓ₅` 和 `Ṁₓ` 的部分表达式
`px0`	`pₓ₀`	多个 `kₓₙ` 表达式和 `Ṁₓ` 的部分表达式
`px1`	`pₓ₁`	多个 `kₓₙ` 表达式和 `İₓ` 的部分表达式
`px2`	`pₓ₂`	多个 `kₓₙ` 表达式和 `ω̇ₓ` 的部分表达式
`px3`	`pₓ₃`	多个 `kₓₙ` 表达式和 `ėₓ` 的部分表达式
`kx0`	`kₓ₀`	`Fₓ₂` 为 `δeₓ` 的系数
`kx1`	`kₓ₁`	`Fₓ₃` 为 `δeₓ` 的系数
`kx2`	`kₓ₂`	`Fₓ₂` 为 `δIₓ` 的系数
`kx3`	`kₓ₃`	`Fₓ₃` 为 `δIₓ` 的系数
`kx4`	`kₓ₄`	`Fₓ₂` 为 `δMₓ` 的系数
`kx5`	`kₓ₅`	`Fₓ₃` 为 `δMₓ` 的系数
`kx6`	`kₓ₆`	`sin fₓ` 为 `δMₓ` 的系数
`kx7`	`kₓ₇`	`Fₓ₂` 为 `pₓ₄` 的系数
`kx8`	`kₓ₈`	`Fₓ₃` 为 `pₓ₄` 的系数
`kx9`	`kₓ₉`	`sin fₓ` 为 `pₓ₄` 的系数
`kx10`	`kₓ₁₀`	`Fₓ₂` 为 `pₓ₅` 的系数
`kx11`	`kₓ₁₁`	`Fₓ₃` 为 `pₓ₅` 的系数
`第三体点.倾角`	`İₓ`	太阳（`İₛ`）或月亮（`İₗ`）对倾角的长期贡献
`第三体赤经点`	`Ω̇ₓ`	太阳（`Ω̇ₛ`）或月亮（`Ω̇ₗ`）对升交点的赤经的长期贡献
`第三体点.eccentricity`	`ėₓ`	太阳（`ėₛ`）或月球（`ėₗ`）对偏心率的长期贡献
`第三体点.近地点幅角`	`ω̇ₓ`	太阳（`ω̇ₛ`）或月球（`ω̇ₗ`）对近地点幅角的长期贡献
`第三体点.mean_anomaly`	`Ṁₓ`	太阳（`Ṁₛ`）或月球（`Ṁₗ`）对平均经度的长期贡献

第三体传播变量

以下变量依赖于传播时间 t。

变量	符号	描述
`第三体平均经度`	`Mₓ`	太阳（`Mₛ`）或月球（`Mₗ`）的平均经度
`fx`	`fₓ`	第三体真近点角
`fx2`	`Fₓ₂`	第三体长周期周期性贡献的部分表达式
`fx3`	`Fₓ₃`	第三体长周期周期性贡献的部分表达式
`第三体偏心率变化`	`δeₓ`	太阳（`δeₛ`）或月球（`δeₗ`）对偏心率的长期周期性贡献
`第三体倾角变化`	`δIₓ`	太阳（`δIₛ`）或月球（`δIₗ`）对倾角的长期周期性贡献
`第三体平均运动变化`	`δMₓ`	太阳（`δMₛ`）或月球（`δMₗ`）对平均运动的长期周期性贡献
`px4`	`pₓ₄`	太阳（`pₛ₄`）或月球（`pₗ₄`）对升交点经度和近地点幅角的长周期周期性贡献的部分表达式
`px5`	`pₓ₅`	太阳（`pₛ₅`）或月球（`pₗ₅`）对升交点经度的长周期周期性贡献的部分表达式

数学表达式

UT1 到儒略历的转换

纪元（自UTC 2000年1月1日12:00以来的儒略年）可以通过AFSPC公式计算

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

或者同一公式的更精确版本

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⁹)

常见初始化

a₁ = (kₑ / n₀)²ᐟ³

      3      3 cos²I₀ - 1
 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

近地初始化

仅在 n₀ > 2π / 225（近地）时定义。

       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⁻¹²

高空近地初始化

仅在 n₀ > 2π / 225（近地）且 p₃ ≥ 220 / (aₑ + 1)（高空）时定义。

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₁²))

椭圆高空近地初始化

仅在 n₀ > 2π / 225（近地）、p₃ ≥ 220 / (aₑ + 1)（高空）且 e₀ > 10⁻⁴（椭圆）时定义。

                    J₃ p₇ ξ  n₀" sin I₀
k₁₂ = - 2 B* cos ω₀ -- ----------------
                    J₂        e₀

        2 p₇ B*
k₁₃ = - - -----
        3 e₀ η

深空初始化

仅在 n₀ ≤ 2π / 225（深空）时定义。

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 - cos²Iₗ)¹ᐟ²

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π

第三体摄动

仅在 n₀ ≤ 2π / 225（深空）时定义。

以下变量为两个第三体（太阳（太阳摄动 s）和月球（月球摄动 l））评估。特定于第三体的变量用 x 注释。在其他部分，x 是 s 或 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₀²)

共振深空初始化

仅在 n₀ ≤ 2π / 225（深空）且以下任一条件满足时定义：

0.0034906585 < n₀ < 0.0052359877（地球同步）
8.26 × 10⁻³ ≤ n₀" ≤ 9.24 × 10⁻³ 和 e₀ ≥ 0.5 （摩尼雅）

恒星时 θ₀ 可以使用国际天文学联合会公式计算

c₂₀₀₀ = y₂₀₀₀ / 100

θ₀ = ¹/₂₄₀ (π / 180) (- 6.2 × 10⁻⁶ c₂₀₀₀³ + 0.093104 c₂₀₀₀²
     + (876600 × 3600 + 8640184.812866) c₂₀₀₀ + 67310.54841) mod 2π

或美国空军航天司令部公式计算

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

地球同步深空初始化

仅在 n₀" ≤ 2π / 225（深空）和 0.0034906585 < n₀" < 0.0052359877（地球同步轨道）时定义。

p₁₇ = 3 (n / a₀")²

𝛿ᵣ₁ = p₁₇ (¹⁵/₁₆ sin²I₀ (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₀"²

摩尼亞深空初始化

仅在 n₀" ≤ 2π / 225（深空）和 8.26 × 10⁻³ ≤ n₀" ≤ 9.24 × 10⁻³ 以及 e₀ ≥ 0.5（摩尼雅）时定义。

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⁻⁶ (³/₂ sin²I₀) 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 sin²I₀ F₂₂₀) G₄₁₀

D₄₄₂₂ = 2 p₂₀ 7.3636953 × 10⁻⁹ (³¹⁵/₈ sin⁴I₀) G₄₂₂

D₅₂₂₀ = p₂₁ 1.1428639 × 10⁻⁷ (³¹⁵/₃₂ sin I₀
        (sin²I₀ (1 - 2 p₁ - 5 p₁²)
        + 0.33333333 (- 2 + 4 p₁ + 6 p₁²))) G₅₂₀

D₅₂₃₂ = p₂₁ 1.1428639 × 10⁻⁷ (sin I₀
        (4.92187512 sin²I₀ (- 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₅₃₃

常见传播

以下值取决于传播时间 t（自原点以来的分钟数）。

命名条件具有以下含义

near earth： n₀" ≤ 2π / 225
low altitude near earth： near earth 和 p₃ < 220 / (aₑ + 1)
high altitude near earth： near earth 和 p₃ ≥ 220 / (aₑ + 1)
elliptic high altitude near earth： high altitude near earth 和 e₀ > 10⁻⁴
non-elliptic near earth： low altitude near earth 或 high altitude near earth 和 e₀ ≤ 10⁻⁴
deep space： n₀" > 2π / 225
non-Lyddane deep space： deep space 和 I ≥ 0.2
Lyddane deep space： deep space 和 I < 0.2
AFSPC Lyddane deep space： Lyddane deep space 并使用与原始AFSPC实现相同的表达式，在 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 - cos²I othewise

p₃₄ = │ k₄          if near earth
      │ 7 cos²I - 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 cos²I - 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

近地传播

仅在 n₀ > 2π / 225（近地）时定义。

p₂₄ = M₀ + Ṁ t

p₂₇ = | e₀ - (C₄ t + C₅ (sin p₂₆ - k₇)) if high altitude
      | e₀ - C₄ t                       otherwise

高空近地传播

仅在 n₀ > 2π / 225（近地）且 p₃ ≥ 220 / (aₑ + 1)（高空）时定义。

elliptic 表示 e₀ > 10⁻⁴。

p₂₅ = k₁₃ ((1 + η cos p₂₄)³ - k₁₁) + k₁₂ t

p₂₆ = │ p₂₄ + p₂₅ if elliptic
      │ p₂₄       otherwise

深空传播

仅在 n₀ ≤ 2π / 225（深空）时定义。

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

第三体传播

仅在 n₀ ≤ 2π / 225（深空）时定义。

以下变量为两个第三体（太阳（太阳摄动 s）和月球（月球摄动 l））评估。特定于第三体的变量用 x 注释。在其他部分，x 是 s 或 l。

Mₓ = Mₓ₀ + nₓ t

fₓ = Mₓ + 2 eₓ sin Mₓ

Fₓ₂ = ¹/₂ sin²fₓ - ¹/₄

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ₓ₃

共振深空传播

仅在 n₀ ≤ 2π / 225（深空）且以下任一条件满足时定义：

0.0034906585 < n₀ < 0.0052359877（地球同步）
8.26 × 10⁻³ ≤ n₀" ≤ 9.24 × 10⁻³ 和 e₀ ≥ 0.5 （摩尼雅）

θ = θ₀ + 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 深空传播

仅在 n₀" ≤ 2π / 225（深空）和 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₂₂

参考文献

大卫·A·瓦拉多、保罗·克劳福德、R·S·胡斯克和T·S·凯尔索，《重访太空轨迹报告#3》，在AIAA/AAS航天动力学专家会议上发表，科罗拉多州科伊塞斯特，2006年8月21日至24日，https://doi.org/10.2514/6.2006-6753

F·R·胡茨、P·W·舒马赫小和R·A·格洛弗，《美国空间监视系统中分析轨道建模的历史》，制导、控制和动力学杂志，第27卷第2期，第174-185页，2004年，https://doi.org/10.2514/1.9161

F·R·胡茨和R·L·罗伊奇，《太空轨迹报告No. 3：NORAD元素集传播模型》，美国空军太空防御司令部，科罗拉多斯普林斯，科罗拉多州，1980年，https://www.celestrak.com/NORAD/documentation/

R·S·胡斯克，《无阻力的共振卫星的有限四体解》，太空轨迹项目报告1，美国空军太空防御司令部，科罗拉多斯普林斯，科罗拉多州，1979年11月，https://doi.org/10.21236/ada081263

依赖关系

~1-1.7MB
~29K SLoC