docs: bilingual comments for the orbit validation scripts

This commit is contained in:
ApsisRS contributors 2026-10-02 01:42:39 +08:00
parent 78614d0a97
commit 1b0d288d73
3 changed files with 156 additions and 13 deletions

View File

@ -33,6 +33,39 @@ Usage:
Exit status is zero when every quantity agrees within tolerance and one Exit status is zero when every quantity agrees within tolerance and one
otherwise. 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 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 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 # 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. # radians/degrees slip, or a dropped rotation, not to allow a loose match.
# 相符程度在浮点水平,因为两侧求的是相同闭式。这些容差是为捕获符号写反、
# 弧度与度混淆,或漏掉的转动而设,并非为允许宽松匹配。
LENGTH_TOLERANCE_KM = 1e-9 LENGTH_TOLERANCE_KM = 1e-9
VELOCITY_TOLERANCE_KMS = 1e-11 VELOCITY_TOLERANCE_KMS = 1e-11
ANGLE_TOLERANCE_RAD = 1e-12 ANGLE_TOLERANCE_RAD = 1e-12
@ -53,11 +88,15 @@ RELATIVE_TOLERANCE = 1e-12
DEFAULT_INPUT = Path("validation/data/expected/rust_kepler.csv") DEFAULT_INPUT = Path("validation/data/expected/rust_kepler.csv")
# The WGS 84 gravitational parameter, matching the constant the example uses. # The WGS 84 gravitational parameter, matching the constant the example uses.
# WGS 84 引力参数,与示例所用的常量一致。
MU = 398_600.5 MU = 398_600.5
def rotation_matrix(inclination: float, right_ascension: float, argument_of_perigee: float): 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) sw, cw = math.sin(argument_of_perigee), math.cos(argument_of_perigee)
si, ci = math.sin(inclination), math.cos(inclination) si, ci = math.sin(inclination), math.cos(inclination)
so, co = math.sin(right_ascension), math.cos(right_ascension) 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): 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)] 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): 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) p = a * (1.0 - e * e)
radius = p / (1.0 + e * math.cos(nu)) radius = p / (1.0 + e * math.cos(nu))
x = radius * 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): 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) return 2.0 * math.pi * math.sqrt(a**3 / mu)
def check_states(rows): 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 = { worst = {
"position": 0.0, "position": 0.0,
"velocity": 0.0, "velocity": 0.0,
@ -123,6 +174,7 @@ def check_states(rows):
# The speed is checked through the vis-viva equation, which is # The speed is checked through the vis-viva equation, which is
# independent of the direction, and the energy through its own # independent of the direction, and the energy through its own
# closed form. # closed form.
# 速率经活力公式核对,后者与方向无关;能量经其自身闭式核对。
speed = math.sqrt(sum(component * component for component in rust_velocity)) 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)) 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) 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) worst["energy"] = max(worst["energy"], energy_error)
_ = speed # the magnitude check is folded into the vis-viva step _ = speed # the magnitude check is folded into the vis-viva step
# 速率大小核对已并入活力公式那一步。
return worst return worst
def check_recovery(rows): 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 = { worst = {
"semi_major_axis": 0.0, "semi_major_axis": 0.0,
"eccentricity": 0.0, "eccentricity": 0.0,
@ -152,6 +208,7 @@ def check_recovery(rows):
} }
for row in rows: for row in rows:
# The input orbit for the second table is fixed by the example. # The input orbit for the second table is fixed by the example.
# 第二张表的输入轨道由示例固定。
expected_a = 7_000.0 expected_a = 7_000.0
expected_e = 0.3 expected_e = 0.3
expected_i = 0.7 expected_i = 0.7
@ -160,6 +217,7 @@ def check_recovery(rows):
nu = float(row["true_anomaly_deg"]) nu = float(row["true_anomaly_deg"])
# The recovered semi-major axis and eccentricity must be unchanged. # The recovered semi-major axis and eccentricity must be unchanged.
# 所恢复的半长轴与偏心率必须不变。
worst["semi_major_axis"] = max( worst["semi_major_axis"] = max(
worst["semi_major_axis"], abs(float(row["recovered_semi_major_axis_km"]) - expected_a) 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 # The angles must be reproduced, modulo a full revolution. The
# argument of perigee and the node are compared through the total # argument of perigee and the node are compared through the total
# longitude, which is the quantity that is actually determined. # longitude, which is the quantity that is actually determined.
# 各角必须在相差整数圈的意义下被复现。近地点幅角与升交点经度经总
# 黄经比较,后者才是实际确定的量。
recovered_lon = float(row["recovered_argument_of_perigee_rad"]) + float( recovered_lon = float(row["recovered_argument_of_perigee_rad"]) + float(
row["recovered_right_ascension_rad"] row["recovered_right_ascension_rad"]
) )
@ -185,7 +245,10 @@ def check_recovery(rows):
def read_blocks(path: Path): 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() text = path.read_text().splitlines()
blocks: list[list[str]] = [[]] blocks: list[list[str]] = [[]]
for line in text: for line in text:

View File

@ -13,6 +13,20 @@ Usage:
Exit status is zero when every state agrees within tolerance and one Exit status is zero when every state agrees within tolerance and one
otherwise. 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 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 # velocity. They are far above the observed agreement, which is at the
# floating point level, and far below the accuracy of the SGP4 theory # 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. # itself, which is on the order of a kilometer per day of propagation.
#
# 交叉核对的容差。
#
# 这些是本工作区对所传播位置与速度承诺的限值。它们远高于所观测到的相符
# 程度(浮点水平),也远低于 SGP4 理论自身的精度(约为每传播一天一千米)。
POSITION_TOLERANCE_KM = 1e-3 POSITION_TOLERANCE_KM = 1e-3
VELOCITY_TOLERANCE_KMS = 1e-6 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 # The element sets, in the same order as the Rust example. They are the
# historical test cases distributed with the reference implementation. # historical test cases distributed with the reference implementation.
# 各根数集,与 Rust 示例相同的次序。它们是随参考实现分发的历史测试用例。
ELEMENT_SETS = [ ELEMENT_SETS = [
( (
"1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927", "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]]: 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: with path.open(newline="") as handle:
reader = csv.DictReader(handle) reader = csv.DictReader(handle)
return [row for row in reader if row.get("set")] return [row for row in reader if row.get("set")]
def reference_state(satellite: Satrec, minutes: float): 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) error, position, velocity = satellite.sgp4_tsince(minutes)
if error != 0: if error != 0:
return None return None
@ -109,6 +135,7 @@ def main(argv: list[str]) -> int:
if rust_error and reference is None: if rust_error and reference is None:
# Both sides reject the instant, which is agreement. # Both sides reject the instant, which is agreement.
# 两侧都拒绝该瞬时,这就是一致。
continue continue
if rust_error or reference is None: if rust_error or reference is None:
mismatched_errors += 1 mismatched_errors += 1

View File

@ -31,6 +31,34 @@ Usage:
pixi run python validation/python/orbit/check_tle.py [path] pixi run python validation/python/orbit/check_tle.py [path]
Exit status is zero when every field agrees and one otherwise. 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 from __future__ import annotations
@ -43,6 +71,9 @@ from pathlib import Path
# agree to the last digit. The tolerance is therefore exact, not # agree to the last digit. The tolerance is therefore exact, not
# approximate, and only guards the exponent fields against a difference in # approximate, and only guards the exponent fields against a difference in
# how the two languages round a very small value. # how the two languages round a very small value.
# 该读取器处理的是解析器所见的相同字符,故二者必须精确到最后一位数字。
# 因此该容差是精确的而非近似的,只是为防护指数字段在两种语言对极小值舍入
# 方式上的差异。
ABSOLUTE_TOLERANCE = 1e-12 ABSOLUTE_TOLERANCE = 1e-12
EXPONENT_TOLERANCE = 1e-18 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 # The element sets, in the order the example writes them. The lines are the
# same ones the example parses. # same ones the example parses.
# 各根数集,按示例写出它们的次序。这些行就是示例所解析的那些。
ELEMENT_SETS = { ELEMENT_SETS = {
"vanguard_1": ( "vanguard_1": (
"1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753", "1 00005U 58002B 00179.78495062 .00000023 00000-0 28098-4 0 4753",
@ -87,7 +119,10 @@ ELEMENT_SETS = {
def checksum(line: str) -> int: 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 total = 0
for character in line[:68]: for character in line[:68]:
if character == "-": if character == "-":
@ -98,7 +133,10 @@ def checksum(line: str) -> int:
def field(line: str, start: int, end: int) -> str: 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] 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 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 ``+`` or ``-`` before the mantissa is the sign of the value, and the sign
before the exponent is the sign of the exponent. before the exponent is the sign of the exponent.
读取 NORAD 的隐含小数点指数表示。
平运动的二阶导数与阻力项在存储时,小数点隐含在第一位数字之前,指数不写
``E``:``-11606-4`` 即 ``-0.11606e-4``。尾数之前的 ``+`` 或 ``-`` 是值
的符号,指数之前的符号是指数的符号。
""" """
text = text.strip() text = text.strip()
if not text: if not text:
return 0.0 return 0.0
# The exponent is introduced by the last sign, once a mantissa digit is # The exponent is introduced by the last sign, once a mantissa digit is
# in hand, so scan from the end. # in hand, so scan from the end.
# 指数由最后一个符号引入(此时已有尾数数字在手),故从末尾扫描。
index = len(text) - 1 index = len(text) - 1
while index > 0 and text[index] not in "+-": while index > 0 and text[index] not in "+-":
index -= 1 index -= 1
@ -137,12 +182,16 @@ def parse_exponential(text: str) -> float:
def parse_line_one(line: str) -> dict: 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() epoch = field(line, 18, 32).strip()
year = int(epoch[:2]) year = int(epoch[:2])
fraction = float(epoch[2:]) fraction = float(epoch[2:])
# The two digit year expands with the catalogue convention: 57 to 99 is # The two digit year expands with the catalogue convention: 57 to 99 is
# the twentieth century and 00 to 56 the twenty first. # the twentieth century and 00 to 56 the twenty first.
# 两位年份按目录约定展开:57 至 99 是二十世纪,00 至 56 是二十一世纪。
century = 1900 if year >= 57 else 2000 century = 1900 if year >= 57 else 2000
return { return {
"satellite_number": int(field(line, 2, 7)), "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: 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 # The eccentricity is stored as a leading decimal: the column holds the
# fractional digits and the decimal point is implied before the first. # fractional digits and the decimal point is implied before the first.
# 偏心率以隐含的前导小数点存储:该列保存小数位,小数点隐含在第一位之前。
eccentricity = float("0." + field(line, 26, 33).strip()) eccentricity = float("0." + field(line, 26, 33).strip())
return { return {
"inclination_deg": float(field(line, 8, 16)), "inclination_deg": float(field(line, 8, 16)),