diff --git a/validation/python/orbit/check_kepler.py b/validation/python/orbit/check_kepler.py index fb929c3..4bc2bd0 100644 --- a/validation/python/orbit/check_kepler.py +++ b/validation/python/orbit/check_kepler.py @@ -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: diff --git a/validation/python/orbit/check_sgp4.py b/validation/python/orbit/check_sgp4.py index f113eb4..76f4c5a 100644 --- a/validation/python/orbit/check_sgp4.py +++ b/validation/python/orbit/check_sgp4.py @@ -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 diff --git a/validation/python/orbit/check_tle.py b/validation/python/orbit/check_tle.py index 8de85e8..0ca1eef 100644 --- a/validation/python/orbit/check_tle.py +++ b/validation/python/orbit/check_tle.py @@ -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)),