ApsisRS
/
apsis-rs
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

docs: bilingual comments for the orbit validation scripts

This commit is contained in:
ApsisRS contributors 2026-10-02 01:42:39 +08:00
parent 8e2a956fa5
commit 0e3d735935
3 changed files with 156 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
validation/python/orbit/check_kepler.py
View File

@ -33,6 +33,39 @@ Usage:
Exit status is zero when every quantity agrees within tolerance and one
otherwise.
用一份独立计算交叉核对 ApsisRS 的二体模型。
读取由 ``cargo run -p apsis-rs --features orbit --example dump_kepler``
写出的量，并只用标准库与 numpy 从轨道力学的闭式表达式重算。要点在于
Python 侧是依据几何而非依据 Rust 代码写成的，故 Rust 公式中的错误会表现
为不一致。
所用的关系都是标准的：
* 真近点角 ``nu`` 处的近焦点坐标位置与速度，
p = a (1 - e^2), r = p / (1 + e cos nu)
x = r cos nu, y = r sin nu
vx = -sqrt(mu / p) sin nu
vy = sqrt(mu / p) (e + cos nu)
由 ``R3(O) R1(i) R3(w)`` 转入惯性系；
* 活力公式 ``v^2 = mu (2/r - 1/a)``，用于独立于方向核对速率；
* 比轨道能量 ``E = v^2/2 - mu/r``；
* 开普勒第三定律给出的周期 ``T = 2 pi sqrt(a^3/mu)``；
* 圆锥曲线半径 ``r = p / (1 + e cos nu)``。
第二张表按相反方向核对：把状态转回根数，所恢复的根数必须等于产生它的
那些。
本脚本不需要网络访问。
用法：
pixi run python validation/python/orbit/check_kepler.py [path]
每个量都在容差内一致时退出状态为零，否则为一。
"""
from __future__ import annotations
@ -45,6 +78,8 @@ from pathlib import Path
# The agreement is at the floating point level because both sides evaluate
# the same closed forms. The tolerances are set to catch a transposed sign, a
# radians/degrees slip, or a dropped rotation, not to allow a loose match.
# 相符程度在浮点水平，因为两侧求的是相同闭式。这些容差是为捕获符号写反、
# 弧度与度混淆，或漏掉的转动而设，并非为允许宽松匹配。
LENGTH_TOLERANCE_KM = 1e-9
VELOCITY_TOLERANCE_KMS = 1e-11
ANGLE_TOLERANCE_RAD = 1e-12
@ -53,11 +88,15 @@ RELATIVE_TOLERANCE = 1e-12
DEFAULT_INPUT = Path("validation/data/expected/rust_kepler.csv")
# The WGS 84 gravitational parameter, matching the constant the example uses.
# WGS 84 引力参数，与示例所用的常量一致。
MU = 398_600.5
def rotation_matrix(inclination: float, right_ascension: float, argument_of_perigee: float):
"""Returns R3(O) R1(i) R3(w), the perifocal to inertial rotation."""
"""Returns R3(O) R1(i) R3(w), the perifocal to inertial rotation.
返回 R3(O) R1(i) R3(w)，即近焦点系到惯性系的转动。
"""
sw, cw = math.sin(argument_of_perigee), math.cos(argument_of_perigee)
si, ci = math.sin(inclination), math.cos(inclination)
so, co = math.sin(right_ascension), math.cos(right_ascension)
@ -69,12 +108,18 @@ def rotation_matrix(inclination: float, right_ascension: float, argument_of_peri
def apply(matrix, vector):
"""Applies a 3x3 matrix to a vector."""
"""Applies a 3x3 matrix to a vector.
把 3x3 矩阵施加于一个矢量。
"""
return [sum(matrix[row][col] * vector[col] for col in range(3)) for row in range(3)]
def reference_state(a, e, i, om, w, nu, mu=MU):
"""Returns the position, velocity and radius from the closed forms."""
"""Returns the position, velocity and radius from the closed forms.
由闭式返回位置、速度与半径。
"""
p = a * (1.0 - e * e)
radius = p / (1.0 + e * math.cos(nu))
x = radius * math.cos(nu)
@ -91,12 +136,18 @@ def reference_state(a, e, i, om, w, nu, mu=MU):
def reference_period(a, mu=MU):
"""Returns the orbital period from Kepler's third law."""
"""Returns the orbital period from Kepler's third law.
由开普勒第三定律返回轨道周期。
"""
return 2.0 * math.pi * math.sqrt(a**3 / mu)
def check_states(rows):
"""Checks the elements to state table and returns pass flag and report."""
"""Checks the elements to state table and returns pass flag and report.
检查「根数到状态」表，并返回通过标志与报告。
"""
worst = {
"position": 0.0,
"velocity": 0.0,
@ -123,6 +174,7 @@ def check_states(rows):
# The speed is checked through the vis-viva equation, which is
# independent of the direction, and the energy through its own
# closed form.
# 速率经活力公式核对，后者与方向无关；能量经其自身闭式核对。
speed = math.sqrt(sum(component * component for component in rust_velocity))
vis_viva = math.sqrt(mu := MU * (2.0 / math.hypot(*rust_position) - 1.0 / a))
energy_reference = vis_viva * vis_viva / 2.0 - MU / math.hypot(*rust_position)
@ -136,12 +188,16 @@ def check_states(rows):
worst["energy"] = max(worst["energy"], energy_error)
_ = speed # the magnitude check is folded into the vis-viva step
# 速率大小核对已并入活力公式那一步。
return worst
def check_recovery(rows):
"""Checks the state to elements table and returns the worst deviations."""
"""Checks the state to elements table and returns the worst deviations.
检查「状态到根数」表，并返回最差偏差。
"""
worst = {
"semi_major_axis": 0.0,
"eccentricity": 0.0,
@ -152,6 +208,7 @@ def check_recovery(rows):
}
for row in rows:
# The input orbit for the second table is fixed by the example.
# 第二张表的输入轨道由示例固定。
expected_a = 7_000.0
expected_e = 0.3
expected_i = 0.7
@ -160,6 +217,7 @@ def check_recovery(rows):
nu = float(row["true_anomaly_deg"])
# The recovered semi-major axis and eccentricity must be unchanged.
# 所恢复的半长轴与偏心率必须不变。
worst["semi_major_axis"] = max(
worst["semi_major_axis"], abs(float(row["recovered_semi_major_axis_km"]) - expected_a)
)
@ -169,6 +227,8 @@ def check_recovery(rows):
# The angles must be reproduced, modulo a full revolution. The
# argument of perigee and the node are compared through the total
# longitude, which is the quantity that is actually determined.
# 各角必须在相差整数圈的意义下被复现。近地点幅角与升交点经度经总
# 黄经比较，后者才是实际确定的量。
recovered_lon = float(row["recovered_argument_of_perigee_rad"]) + float(
row["recovered_right_ascension_rad"]
)
@ -185,7 +245,10 @@ def check_recovery(rows):
def read_blocks(path: Path):
"""Returns the two CSV blocks separated by a blank line."""
"""Returns the two CSV blocks separated by a blank line.
返回以空行分隔的两个 CSV 块。
"""
text = path.read_text().splitlines()
blocks: list[list[str]] = [[]]
for line in text:

31
validation/python/orbit/check_sgp4.py
View File

@ -13,6 +13,20 @@ Usage:
Exit status is zero when every state agrees within tolerance and one
otherwise.
用参考库交叉核对 ApsisRS 的 SGP4 传播器。
读取由 ``cargo run -p apsis-rs --features orbit --example dump_orbit`` 写出的
状态，用 ``sgp4`` 参考实现传播同一批根数集，并报告每个根数集的最大位置与
速度差异。
本脚本不需要网络访问。
用法：
pixi run python validation/python/orbit/check_sgp4.py [path]
每个状态都在容差内一致时退出状态为零，否则为一。
"""
from __future__ import annotations
@ -30,6 +44,11 @@ from sgp4.model import Satrec, WGS72
# velocity. They are far above the observed agreement, which is at the
# floating point level, and far below the accuracy of the SGP4 theory
# itself, which is on the order of a kilometer per day of propagation.
#
# 交叉核对的容差。
#
# 这些是本工作区对所传播位置与速度承诺的限值。它们远高于所观测到的相符
# 程度（浮点水平），也远低于 SGP4 理论自身的精度（约为每传播一天一千米）。
POSITION_TOLERANCE_KM = 1e-3
VELOCITY_TOLERANCE_KMS = 1e-6
@ -37,6 +56,7 @@ DEFAULT_INPUT = Path("validation/data/expected/rust_sgp4.csv")
# The element sets, in the same order as the Rust example. They are the
# historical test cases distributed with the reference implementation.
# 各根数集，与 Rust 示例相同的次序。它们是随参考实现分发的历史测试用例。
ELEMENT_SETS = [
(
"1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927",
@ -58,14 +78,20 @@ ELEMENT_SETS = [
def read_rows(path: Path) -> list[dict[str, str]]:
"""Reads the Rust output rows."""
"""Reads the Rust output rows.
读取 Rust 输出各行。
"""
with path.open(newline="") as handle:
reader = csv.DictReader(handle)
return [row for row in reader if row.get("set")]
def reference_state(satellite: Satrec, minutes: float):
"""Returns the reference position and velocity, or None on error."""
"""Returns the reference position and velocity, or None on error.
返回参考位置与速度，出错时返回 None。
"""
error, position, velocity = satellite.sgp4_tsince(minutes)
if error != 0:
return None
@ -109,6 +135,7 @@ def main(argv: list[str]) -> int:
if rust_error and reference is None:
# Both sides reject the instant, which is agreement.
# 两侧都拒绝该瞬时，这就是一致。
continue
if rust_error or reference is None:
mismatched_errors += 1

61
validation/python/orbit/check_tle.py
View File

@ -31,6 +31,34 @@ Usage:
pixi run python validation/python/orbit/check_tle.py [path]
Exit status is zero when every field agrees and one otherwise.
用一份独立的读取器交叉核对 ApsisRS 的两行根数解析器。
读取由 ``cargo run -p apsis-rs --features orbit --example dump_tle`` 写出的
根数字段，并用一个依据 NORAD 格式写成的按列索引读取器重新解析同一批根数
集，然后比较每个字段。每一行的校验和都从其定义重算，故字段放错列或符号
约定被丢弃会表现为不一致。
该格式是按位置的，行的第一个字符在零列：
* 第 1 行：卫星编号在第 2 至 6 列，密级在第 7 列，国际编号在第 9 至 16
列，历元年与一天中的小数部分在第 18 至 32 列，平运动的一阶导数除以二
在第 33 至 43 列，二阶导数除以六在第 44 至 52 列，阻力项在第 53 至 61
列，根数集编号在第 62 至 67 列，校验和在第 68 列；
* 第 2 行：倾角在第 8 至 16 列，赤经在第 17 至 25 列，偏心率在第 26 至 33
列以隐含的前导小数点表示，近地点幅角在第 34 至 42 列，平近点角在第 43
至 51 列，平运动在第 52 至 63 列，圈数在第 63 至 68 列，校验和在第 68
列。
连字符计为一，每个数字计为其自身值，在前六十八列上求和并对十取模。
本脚本不需要网络访问。
用法：
pixi run python validation/python/orbit/check_tle.py [path]
每个字段都一致时退出状态为零，否则为一。
"""
from __future__ import annotations
@ -43,6 +71,9 @@ from pathlib import Path
# agree to the last digit. The tolerance is therefore exact, not
# approximate, and only guards the exponent fields against a difference in
# how the two languages round a very small value.
# 该读取器处理的是解析器所见的相同字符，故二者必须精确到最后一位数字。
# 因此该容差是精确的而非近似的，只是为防护指数字段在两种语言对极小值舍入
# 方式上的差异。
ABSOLUTE_TOLERANCE = 1e-12
EXPONENT_TOLERANCE = 1e-18
@ -50,6 +81,7 @@ DEFAULT_INPUT = Path("validation/data/expected/rust_tle.csv")
# The element sets, in the order the example writes them. The lines are the
# same ones the example parses.
# 各根数集，按示例写出它们的次序。这些行就是示例所解析的那些。
ELEMENT_SETS = {
"vanguard_1": (
"1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753",
@ -87,7 +119,10 @@ ELEMENT_SETS = {
def checksum(line: str) -> int:
"""Returns the modulo ten checksum of a line from its definition."""
"""Returns the modulo ten checksum of a line from its definition.
由其定义返回某行对十取模的校验和。
"""
total = 0
for character in line[:68]:
if character == "-":
@ -98,7 +133,10 @@ def checksum(line: str) -> int:
def field(line: str, start: int, end: int) -> str:
"""Returns a slice of a line."""
"""Returns a slice of a line.
返回某行的一个切片。
"""
return line[start:end]
@ -110,12 +148,19 @@ def parse_exponential(text: str) -> float:
written without an ``E``: ``-11606-4`` is ``-0.11606e-4``. A leading
``+`` or ``-`` before the mantissa is the sign of the value, and the sign
before the exponent is the sign of the exponent.
读取 NORAD 的隐含小数点指数表示。
平运动的二阶导数与阻力项在存储时，小数点隐含在第一位数字之前，指数不写
``E``：``-11606-4`` 即 ``-0.11606e-4``。尾数之前的 ``+`` 或 ``-`` 是值
的符号，指数之前的符号是指数的符号。
"""
text = text.strip()
if not text:
return 0.0
# The exponent is introduced by the last sign, once a mantissa digit is
# in hand, so scan from the end.
# 指数由最后一个符号引入（此时已有尾数数字在手），故从末尾扫描。
index = len(text) - 1
while index > 0 and text[index] not in "+-":
index -= 1
@ -137,12 +182,16 @@ def parse_exponential(text: str) -> float:
def parse_line_one(line: str) -> dict:
"""Parses the first line by column index."""
"""Parses the first line by column index.
按列索引解析第一行。
"""
epoch = field(line, 18, 32).strip()
year = int(epoch[:2])
fraction = float(epoch[2:])
# The two digit year expands with the catalogue convention: 57 to 99 is
# the twentieth century and 00 to 56 the twenty first.
# 两位年份按目录约定展开：57 至 99 是二十世纪，00 至 56 是二十一世纪。
century = 1900 if year >= 57 else 2000
return {
"satellite_number": int(field(line, 2, 7)),
@ -159,9 +208,13 @@ def parse_line_one(line: str) -> dict:
def parse_line_two(line: str) -> dict:
"""Parses the second line by column index."""
"""Parses the second line by column index.
按列索引解析第二行。
"""
# The eccentricity is stored as a leading decimal: the column holds the
# fractional digits and the decimal point is implied before the first.
# 偏心率以隐含的前导小数点存储：该列保存小数位，小数点隐含在第一位之前。
eccentricity = float("0." + field(line, 26, 33).strip())
return {
"inclination_deg": float(field(line, 8, 16)),