力模型

e2m2e 的力模型子包提供可配置、可组合的航天器摄动力模型，支持从配置字典构建、序列化到 JSON，并与星历动力学系统配合完成轨道传播。

核心概念

这三个概念构成力模型子包的核心接口：

• PhysicalModel：所有力模型的抽象基类，定义 compute_acceleration(t, state, system) 接口。另有两个可选钩子：compute_jacobian(t, state, system) 返回解析雅可比 ∂a/∂r（默认返回 None，由 ForceModel 有限差分兜底）；to_rust_spec(system) 把力序列化为 Rust 编译路径接受的元组（默认返回 None，表示该力不支持 Rust 编译）。

• ForceModel：力模型组合，把多个 PhysicalModel 组合成一次传播所需的运动方程；按名登记、启用/禁用，并通过 Rust 积分器完成传播。

• ForceEntry：力模型组合内单个力模型的登记记录，包含 name、force 和 enabled 三个字段。

支持力模型类型

类型	说明	配置 type 名
PointMassGravity	中心天体点质量引力	PointMassGravity
ThirdBodyGravity	第三体引力摄动（含间接项）	ThirdBodyGravity
IndirectTerm	单独的间接项修正	IndirectTerm
GravityField	球谐重力场（天体无关，支持 EGM96/GRGM900C 及自定义 .gfc/.cof）	GravityField
DragModel	大气阻力（依赖注入大气密度模型）	DragModel
SolarRadiationPressure	太阳光压（cannonball + 可选阴影模型）	SolarRadiationPressure
EcomSolarRadiationPressure	ECOM 经验光压（9 系数 DYB 参数化，DFH 兼容）	EcomSolarRadiationPressure
FiniteBurn	连续推力（封闭 DSL：constant/pulse 推力 + fixed 方向）	FiniteBurn
VariableMassFiniteBurn	可变质量连续推力（质量作为状态量随推力消耗，7D 受控动力学）	—（不走配置注册表，直接构造）
RelativisticCorrection	后牛顿相对论修正（Schwarzschild / Lense-Thirring / de Sitter）	RelativisticCorrection

另有机动事件 ImpulsiveBurn``（瞬时 Δv，在指定 epoch 直接修改状态速度）， 它不是 ``PhysicalModel，不参与加速度叠加，也不进配置注册表；施加方式见 下文「传播接口」小节的 propagate_maneuvers。

警告

用 GravityField 模拟月球（含 degree=0 中心项）时，必须单独补月球间接项 IndirectTerm("MOON")，且不能再加 ThirdBodyGravity("MOON")——后者会与 GravityField 的中心项重复计算月球点质量。

内置力模型公式

PointMassGravity — 中心天体点质量引力。适用于参考系原点天体自身的二体引力：

\[\mathbf{a} = -\frac{\mu}{|\mathbf{r}|^3} \, \mathbf{r}\]

其中 μ 为引力参数（km³/s²），r 为航天器相对中心天体的位置。 它不查任何第三体的位置，只能表达参考系原点天体自身的引力；其他天体的 引力贡献用 ThirdBodyGravity。

ThirdBodyGravity — 参考系原点之外的天体引力摄动。一个实例对应一个摄动 天体（如 ThirdBodyGravity("MOON")）。加速度由直接项与间接项合成：

\[\mathbf{a} = -\mu_i \left[ \frac{\mathbf{r} - \mathbf{r}_i}{|\mathbf{r} - \mathbf{r}_i|^3} + \frac{\mathbf{r}_i}{|\mathbf{r}_i|^3} \right]\]

其中 r 为航天器相对原点的位置，r_i 为摄动天体相对原点的位置（由 SPICE 查询）。间接项扣除摄动天体对原点的引力，保持坐标原点固定。

GravityField — 完全正规化球谐重力场（Cnm/Snm），用 Pines 递推计算非球形 引力加速度。天体无关：地球（EGM96）、月球（GRGM900C）等共用同一个类，按 body 参数自动切换 body-fixed 轴与系数文件。位势展开为：

\[U = \frac{\mu}{r} \sum_{n=0}^{N} \left(\frac{R}{r}\right)^n \sum_{m=0}^{n} \left(C_{nm}\cos m\lambda + S_{nm}\sin m\lambda\right) \bar{P}_{nm}(\sin\phi)\]

加速度由 Pines 方法直接递推位势梯度得到，不经过球谐系数的解析微分。内置 固体潮修正：地球支持 Step1（天体无关）+ Step2（频率相关）+ 极潮 + 永久潮； 月球支持 k₂ = 0.024116 Love 数的固体潮。

DragModel — 大气阻力。在 ITRF（地固系）中计算密度与相对速度，自动完成 参考系↔ITRF 坐标变换。大气在 ITRF 中静止，相对速度等于航天器 ITRF 速度：

\[\mathbf{a}_{\text{drag}} = -\frac{1}{2} \, \rho \, \frac{C_d A}{m} \, |\mathbf{v}_{\text{rel}}| \, \mathbf{v}_{\text{rel}}\]

其中 ρ 为大气密度（由 ExponentialAtmosphere 提供，US Standard Atmosphere 1976 分段指数模型），C_d 为阻力系数（默认 2.2），A/m 为 面积质量比。

FiniteBurn — 连续推力加速度力模型。推力大小（thrust_profile(t) → 标量 N）与方向（direction）解耦，方向支持传播惯性系、VNB、LVLH 三种 坐标系：

\[\mathbf{a}_{\text{thrust}} = \frac{T(t)}{m} \, \hat{\mathbf{d}}\]

其中 T(t) 为标量推力函数，d̂ 为归一化方向向量。配置往返支持固定 推力/脉冲剖面与固定方向的封闭 DSL；只有该 DSL 控制可在 Rust 编译传播中执行， 任意 Python callable 会在传播入口显式拒绝。

VariableMassFiniteBurn — 可变质量连续推力，是低推力转移与月面动力下降等 最优控制问题的受控动力学基座。与 FiniteBurn 的唯一区别：质量不是常量， 而是状态量 state[6]，随推力按 ṁ = −T/(Isp·g₀) 消耗：

\[\mathbf{a}_{\text{thrust}} = \frac{T}{m} \, \hat{\mathbf{d}}, \qquad \dot{m} = -\frac{T}{I_{\text{sp}} \, g_0}\]

含它的 ForceModel.propagate 把状态扩展为 7D [r, v, m]，自动分流到 Rust propagate_compiled_lowthrust``（受控 EOM 复用 ``augmented_state::augmented_eom_7d）。本期仅支持常量推力与固定方向 （to_rust_spec 非 None 才走 Rust 路径）；可调用方向暂返回 NotImplementedError。

import numpy as np
from e2m2e.algorithm.forces import ForceModel, GravityField, VariableMassFiniteBurn

burn = VariableMassFiniteBurn(
    thrust=0.1,            # N
    isp=3000.0,            # s
    initial_mass=1000.0,   # kg，写入状态第 7 维
    direction=np.array([0.0, 1.0, 0.0]),  # 固定方向（圆轨道初速方向）
)
fm = ForceModel(system, forces=[GravityField("EARTH", degree=0, order=0), burn])

# 7D 初值：[r, v, m]
state0 = np.array([6678.137, 0.0, 0.0, 0.0, 7.726, 0.0, 1000.0])
result = fm.propagate(state0, (et0, et0 + 86400.0))
# result["states"] 形状 (n, 7)，最后一列为质量剖面

RelativisticCorrection — 后牛顿相对论修正，含三项，公式与 GMAT 对齐：

• Schwarzschild 项（质量引起的时空弯曲）：

  \[\mathbf{a}_S = \frac{\gamma \mu}{c^2 r^3} \left[ \left(\frac{4\mu}{r} - v^2\right) \mathbf{r} + 4(\mathbf{r} \cdot \mathbf{v})\mathbf{v} \right]\]

• Lense-Thirring 项（参考系拖曳）：

  \[\mathbf{a}_{LT} = \frac{2\mu}{c^2 r^3} \left[ \frac{3}{r^2}(\mathbf{r} \cdot \mathbf{J})(\mathbf{r} \times \mathbf{v}) + \mathbf{v} \times \mathbf{J} \right]\]

• de Sitter 项（测地进动）：

  \[\mathbf{a}_{dS} = 2 \, \boldsymbol{\omega} \times \mathbf{v}\]

其中 γ = 1.0 为后牛顿参数，c 为光速，J 为天体角动量参数， ω 为 de Sitter 进动角速度。

配置 Schema

顶层配置字典结构：

{
    "version": 1,
    "forces": [
        {
            "name": "j2",
            "type": "GravityField",
            "enabled": True,
            "params": {
                "body": "EARTH",
                "degree": 2,
                "order": 0,
            },
        },
        {
            "name": "drag",
            "type": "DragModel",
            "enabled": True,
            "params": {
                "body": "EARTH",
                "cd": 2.2,
                "area": 10.0,
                "mass": 1000.0,
                "atmosphere": {
                    "type": "ExponentialAtmosphere",
                    "params": {"f107": 150.0, "ap": 4.0},
                },
            },
        },
        {
            "name": "srp",
            "type": "SolarRadiationPressure",
            "enabled": True,
            "params": {
                "area": 5.0,
                "mass": 1000.0,
                "cr": 1.5,
                "shadow": None,  # 全光照
            },
        },
        {
            "name": "burn",
            "type": "FiniteBurn",
            "enabled": False,
            "params": {
                "mass": 1000.0,
                "thrust_profile": {
                    "kind": "constant",
                    "thrust": 0.5,
                },
                "direction": {
                    "kind": "fixed",
                    "vector": [1.0, 0.0, 0.0],
                },
            },
        },
    ],
}

关键字段说明：

• version：当前固定为 1，用于日后 schema 迁移。

• forces：力模型条目数组，顺序即传播时的叠加顺序。

• name：容器内唯一标识，用于 get_force / enable / disable / remove_force。

• type：Python 类名，也是 force_config 注册表的 key。

• enabled：布尔开关；False 时传播跳过该力，但保留在容器内。

• params：构造参数。嵌套依赖（如 atmosphere、shadow）用 {"type": ..., "params": ...} 递归表达；null 表示未注入。

name-based 注册机制

from e2m2e.algorithm.forces import ForceModel, GravityField, DragModel
from e2m2e.algorithm.forces.atmosphere import ExponentialAtmosphere

fm = ForceModel(system)

# 省略 name 时自动按类名生成，遇同类自动消歧
fm.add_force(GravityField("EARTH", degree=2, order=0))
fm.add_force(GravityField("EARTH", degree=4, order=4))   # 自动命名为 GravityField_2

# 显式命名
fm.add_force(DragModel(atmosphere=ExponentialAtmosphere(), area=10.0, mass=1000.0), name="drag")

# 按名操作
fm.disable("drag")          # 暂时关闭阻力
fm.enable("drag")         # 重新打开
fm.get_force("drag")      # 返回 DragModel 实例
fm.remove_force("drag")   # 从容器中移除

# 列出所有注册记录
for entry in fm.list_forces():
    print(f"{entry.name}: {type(entry.force).__name__}, enabled={entry.enabled}")

传播接口

ForceModel.propagate 用 Rust rk_step 单步步进器做自适应传播：

from e2m2e.integrators import RkMethod

result = fm.propagate(
    state0,
    (t0, tf),
    t_eval=t_eval,             # 输出采样点，默认 linspace(t0, tf, 100)
    with_stm=True,             # 同时积分状态转移矩阵
    initial_step=1.0,          # 初始步长，默认从初始状态估算
    events=[...],              # 终止事件列表，见下
    max_steps=200_000,
    method=RkMethod.PD45,      # RK 方法，默认 PD45
)

with_stm=True 时把 6 维状态与 6×6 STM 展平拼成 42 维增广状态一起积分： 各力的解析雅可比 compute_jacobian 直接叠加，不提供解析雅可比的力由 ForceModel 用三点中心差分兜底；STM 分量不参与步长误差控制，接受/拒绝 只看前 6 维物理状态（对齐 GMAT）。返回字典在 time、states、 terminal_event_index 之外多一个 stm 键，形状 (n_points, 6, 6)。

events 是终止事件列表，每个事件是 f(t, state) -> float 的可调用， 函数值符号变化即在该步停止，terminal_event_index 记录触发的事件下标； with_stm=True 时事件函数接收 6 维物理状态（非增广状态）。

带脉冲机动时用 propagate_maneuvers(initial_state, t_span, burns)：按 epoch 排序 ImpulsiveBurn 列表，逐段 coast 传播，在每个 burn epoch 处施加 state[3:6] += delta_v 后续传；返回字典额外含 burns 键， 记录每次机动的施加位置与前后速度。

from e2m2e.algorithm.forces import ImpulsiveBurn

burns = [ImpulsiveBurn(epoch=et0 + 1800.0, delta_v=np.array([0.0, 0.01, 0.0]))]
result = fm.propagate_maneuvers(state0, (et0, et0 + 7200.0), burns)

Rust 编译快速路径

spice feature 启用、无 events、不带 STM，且所有启用力模型的 to_rust_spec() 都非 None 时，propagate 自动分流到 Rust propagate_compiled：力模型一次序列化后整个积分循环在 Rust 内完成， 消除逐步 Python↔Rust 跨界；30 天 NRHO 传播约 9.6 s，Python 路径约 95 s。 任一条件不满足时自动回退 Python 路径，返回格式一致，无需用户干预。

带 STM 时有对应的 Rust compiled STM 快速路径（propagate_compiled_stm_py）， 但 SolarRadiationPressure 与 RelativisticCorrection 无解析雅可比， 含这两个力的 STM 传播不走该路径，回退 Python 增广积分。

Rust 星历预采样缓存

Rust 积分内循环里，GravityField/ThirdBodyGravity/IndirectTerm 每个 RK 子步都跨界调 cspice FFI（spkezr/pxform），即便 GravityField(degree=0) 也不例外。Python 侧的 EphemCache 只拦 Python 层查询，对 Rust 积分内循环无效——Rust 直接走 spk_accel/gravity_field → spice_ffi，不回 Python。

enable_ephem_cache 在积分前把要用到的天体状态与帧旋转矩阵在均匀网格上 预采样、建三次样条存内存表，此后上述力模型每步查表替代 FFI：

from e2m2e._integrators import enable_ephem_cache, disable_ephem_cache

enable_ephem_cache(
    targets=[("MOON", "EARTH"), ("SUN", "EARTH"), ("EARTH", "SOLAR SYSTEM BARYCENTER")],
    frame_pairs=[("ITRF93", "J2000"), ("MOON_PA", "J2000")],
    et_start=et0, et_end=et0 + duration, dt=600.0,
    # 可选：6×6 状态变换对（Lense-Thirring 相对论力需要 body-fixed→J2000）
    sxform_pairs=[("ITRF93", "J2000")],
)
try:
    result = fm.propagate(...)
finally:
    disable_ephem_cache()

三次样条保 C² 连续，避免线性插值网格点导数跳变导致自适应积分器缩步长。 精度与网格步长相关（pxform 600s 网格下末态精度 ~1e-3 km）；未激活时所有 路径走原 FFI，逐字一致。附带的零误差优化：GravityField 在 body == propagation_origin``（地心系地球重力场等常见场景）时跳过 origin→SSB 查询（该量在 ``r_body_icrf 短路分支未被使用）。

配置驱动构建流程

以下示例展示从配置字典构建 ForceModel，接入 EphemerisDynamics 完成一次 LEO 轨道传播。

import numpy as np

from e2m2e.algorithm.coordinate import CoordinateSystem, ICRSAxes
from e2m2e.algorithm.dynamics import EphemerisSystem
from e2m2e.data.kernels.manager import SPICEManager
from e2m2e.algorithm.coordinate.standard_origins import CelestialBodyOrigin
from e2m2e.algorithm.forces import ForceModel

# 1. 准备星历系统（ICRF + 地球中心）
spice = SPICEManager()
spice.load_kernel("path/to/de440.bsp")

system = EphemerisSystem(bodies=["EARTH"], spice=spice, origin="EARTH")
system.coordinate_system = CoordinateSystem(
    axes=ICRSAxes(),
    origin=CelestialBodyOrigin(body="EARTH", spice=spice),
)

# 2. 定义配置字典
config = {
    "version": 1,
    "forces": [
        {
            "name": "j2",
            "type": "GravityField",
            "enabled": True,
            "params": {"body": "EARTH", "degree": 2, "order": 0},
        },
        {
            "name": "drag",
            "type": "DragModel",
            "enabled": True,
            "params": {
                "body": "EARTH",
                "cd": 2.2,
                "area": 10.0,
                "mass": 1000.0,
                "atmosphere": {
                    "type": "ExponentialAtmosphere",
                    "params": {"f107": 150.0, "ap": 4.0},
                },
            },
        },
    ],
}

# 3. 从配置构建 ForceModel
fm = ForceModel.from_config(config, system)

# 4. 设置初始状态（400 km 圆轨道，km / km/s）
r = 6378.137 + 400.0
v = np.sqrt(398600.4415 / r)
state0 = np.array([r, 0.0, 0.0, 0.0, v, 0.0])

# 5. 传播 10 分钟
et0 = spice.utc_to_et("2025-06-21T11:00:06")
t_span = (et0, et0 + 600.0)
t_eval = np.linspace(et0, et0 + 600.0, 50)

result = fm.propagate(state0, t_span, t_eval=t_eval, max_steps=200_000)

print(f"states shape: {result['states'].shape}")   # (50, 6)
print(f"time span: {result['time'][0]} -> {result['time'][-1]}")

# 6. 序列化回配置，再构建一次（验证 round-trip 契约）
# 注意：to_config 输出是规范形式（如 GravityField 恒带 input_frame/gravity_file
# 两键），与手写 config 不一定逐键相等；契约是“再序列化一次字典相等”。
config_roundtrip = fm.to_config()
fm2 = ForceModel.from_config(config_roundtrip, system)
assert fm2.to_config() == config_roundtrip

# 7. 存盘 / 读回
from e2m2e.algorithm.forces import dump_force_config, load_force_config

dump_force_config(fm, "leo_forces.json")
fm_loaded = load_force_config("leo_forces.json", system)

JSON 文件 IO

from e2m2e.algorithm.forces import dump_force_config, load_force_config

# 写入 JSON
dump_force_config(fm, "forces.json")

# 从 JSON 读回
fm2 = load_force_config("forces.json", system)

常见错误与排错

未注册名称

fm.get_force("nonexistent")   # KeyError: 'nonexistent'
fm.disable("nonexistent")      # KeyError: 'nonexistent'

显式命名冲突

fm.add_force(GravityField("EARTH"), name="primary")
fm.add_force(GravityField("EARTH"), name="primary")  # ValueError: force name 'primary' already exists

配置版本不匹配

ForceModel.from_config({"version": 2, "forces": []}, system)
# ValueError: unsupported config version 2; expected 1

未知力类型

from e2m2e.algorithm.forces.force_config import build_force

build_force("UnknownType", {})
# ValueError: unknown force type 'UnknownType'; known types: ['DragModel', 'FiniteBurn', 'GravityField', 'IndirectTerm', 'PointMassGravity', 'RelativisticCorrection', 'SolarRadiationPressure', 'ThirdBodyGravity']

不可序列化的推力剖面

FiniteBurn 若使用用户手写的 lambda 作为 thrust_profile，仍可正常传播，但 to_config() 会抛出 NotSerializableError。解决方式：改用 force_config 提供的 DSL（{"kind": "constant"} 或 {"kind": "pulse"}）构造推力剖面。

力模型容器与 Rust 积分器传播实现。

class e2m2e.algorithm.forces.force_model.ForceEntry(name, force, enabled=True)

基类：object

容器内单个力模型的注册记录。

参数:

• name (str)

• force (PhysicalModel)

• enabled (bool)

name: str

force: PhysicalModel

enabled: bool = True

class e2m2e.algorithm.forces.force_model.ForceModel(system, forces=None)

基类：object

聚合多个 PhysicalModel 并完成传播的动力学容器。

不继承 Dynamics：Dynamics 是 CR3BP/Ephemeris 的基类，其 propagate() 基于 scipy.solve_ivp 与 STM 模板方法；ForceModel 用 Rust rk_step 单步步进器实现自适应传播，支持 with_stm=True （各力解析雅可比叠加、无雅可比的力用有限差分兜底），不支持 Jacobi。 此前形式上继承 Dynamics 只为复用几个数据属性，却全部重写 propagate 并对 STM/Jacobi 抛 NotImplementedError——是 LSP 违反（假继承）。

参数:

• system (Any)

• forces (list[PhysicalModel] | None)

DEFAULT_TOLERANCE = 1e-12

DEFAULT_MAX_STEP = 60.0

STATE_DIM = 6

STM_DIM = 42

__init__(system, forces=None)

初始化 ForceModel。

参数:

• system (Any) -- 动力学系统，必须提供 coordinate_system。

• forces (list[PhysicalModel] | None) -- 初始力模型列表，默认空列表。

返回类型:

None

rtol: float

atol: float

max_step: float

last_trajectory: tuple[ndarray, ndarray] | None

property forces: tuple[PhysicalModel, ...]

当前聚合的力模型，只读（含已 disable 的项）。

add_force(force, name=None)

添加一个力模型。

参数:

• force (PhysicalModel) -- 待添加的力模型。

• name (str | None) -- 力模型的名字。缺省时按类名自动生成，遇同类自动消歧 （Foo、Foo_2、Foo_3…）。显式给出且与已有名字 冲突时抛 ValueError。

返回类型:

None

get_force(name)

按名取力模型；不存在抛 KeyError。

参数:

name (str)

返回类型:

PhysicalModel

list_forces()

返回所有力模型的注册记录（含已 disable 的项）。

与 forces 属性的区别：本方法暴露 name 与 enabled 两个维度。

返回类型:

list[ForceEntry]

enable(name)

按名启用一个力模型；不存在抛 KeyError。

参数:

name (str)

返回类型:

None

disable(name)

按名禁用一个力模型（跳过加速度计算，但保留在容器内）。

不存在抛 KeyError。

参数:

name (str)

返回类型:

None

to_config()

序列化为配置字典 {version, forces: [...]}。

每条力经 force_config.serialize_force 转 {type, params}， 容器补 name 与 enabled。round-trip 契约见 ADR 0004。

返回类型:

dict[str, Any]

classmethod from_config(config, system)

从配置字典构建 ForceModel。

校验 version，逐条 force_config.build_force 构造并按 name 注册；enabled: false 的条目构造后立即 disable。

参数:

• config (dict[str, Any])

• system (Any)

返回类型:

ForceModel

remove_force(index)

移除一个力模型（按索引或名字）。

参数:

index (int | str) -- 整数索引，或力模型注册名。

抛出:

• ValueError -- 名字不存在时。

• IndexError -- 索引越界时。

返回类型:

None

propagate(initial_state, t_span, t_eval=None, with_stm=False, with_jacobi=False, *, initial_step=None, events=None, max_steps=100000, method=None)

使用 Rust 编译传播轨迹（零跨界）。

issue #378：默认传播一律走编译 Rust（propagate_compiled / propagate_compiled_stm_py / propagate_compiled_lowthrust）。 扩展不可用（RustExtensionUnavailableError）或力模型无 Rust spec （NotImplementedError 能力错误）时显式报错，不再静默回退 Python/scipy。

参数:

• initial_state (ArrayLike) -- 初始状态向量，形状 (6,)。

• t_span (tuple[float, float]) -- 时间区间 [t0, tf]，单位为 SPICE et 秒。

• t_eval (ArrayLike | None) -- 评估时间点数组，默认 linspace(t0, tf, 100)。

• with_stm (bool) -- 是否同时积分状态转移矩阵。返回字典额外含 stm 键， 形状 (n_points, 6, 6)。STM 不参与步长误差控制（对齐 GMAT）。

• with_jacobi (bool) -- 不支持，传 True 抛 NotImplementedError。

• initial_step (float | None) -- 初始步长，默认从初始状态估算。

• events (list[Callable[[float, ndarray[tuple[Any, ...], dtype[floating]]], float]] | None) -- 不支持。ForceModel 事件传播需要 compiled-forces Rust API， 当前未提供，传 events 抛 NotImplementedError（不能回退 Python RHS，issue #378）。

• max_steps (int) -- 最大积分步数，默认 100_000。

• method (RkMethod | None) -- Runge-Kutta 积分器方法，默认 PD45。

返回:

包含 time、states 和 terminal_event_index 的字典； with_stm=True 时额外含 stm 键。

返回类型:

dict[str, Any]

propagate_maneuvers(initial_state, t_span, burns, *, initial_step=None, max_steps=100000, method=None)

带脉冲机动的传播：coast 段之间在 burn epoch 处施加 Δv。

按 epoch 排序 burns，依次 coast → 施加 Δv → 续传。burn epoch 处 输出行携带 post-burn 速度（丢 pre-burn 行，无重复 epoch）。

参数:

• initial_state (ArrayLike)

• t_span (tuple[float, float])

• burns (list[ImpulsiveBurn])

• initial_step (float | None)

• max_steps (int)

• method (RkMethod | None)

返回类型:

dict[str, Any]

ForceModel 配置驱动：力模型 ↔ dict 序列化与 JSON IO。

设计见 ADR 0004。容器级编排（信封、version、entry 拼装）在 ForceModel.to_config / from_config；本模块只负责"单条力"的 类型分发与 JSON 文件读写。

e2m2e.algorithm.forces.force_config.serialize_force(force)

把单条力序列化为 {type, params}；未知类型抛 NotSerializableError。

参数:

force (PhysicalModel)

返回类型:

dict[str, Any]

e2m2e.algorithm.forces.force_config.build_force(type_name, params)

按 type 名与 params 构造单条力；未知 type 抛 ValueError。

参数:

• type_name (str)

• params (dict[str, Any])

返回类型:

PhysicalModel

e2m2e.algorithm.forces.force_config.dump_force_config(fm, path)

把 ForceModel.to_config() 的结果写入 JSON 文件。

参数:

• fm (Any)

• path (str | Path)

返回类型:

None

e2m2e.algorithm.forces.force_config.load_force_config(path, system)

从 JSON 文件读取配置并构建 ForceModel。

参数:

• path (str | Path)

• system (Any)

返回类型:

Any

力模型抽象基类。

class e2m2e.algorithm.forces.physical_model.PhysicalModel

基类：object

物理力模型基类。

力模型在 Python 侧只承担"配置定义"职责：参数验证、to_rust_spec 序列化、to_config/from_config。加速度与雅可比计算全部由 Rust 编译路径（ForceModel.propagate → propagate_compiled/ propagate_compiled_stm_py）承载，不保留 Python 参考实现（issue #378）： 需要 Rust 的场景扩展不可用即显式报错，不静默回退到 Python。

所有坐标约定都在 system.coordinate_system 下完成；需要非默认坐标系 计算的子类应通过 system.coordinate_system.transform_state() / transform_vector() 自行完成转换。

to_rust_spec(system)

序列化该 force 为 Rust propagate_compiled 接受的元组。

返回 None 表示该 force 不支持 Rust 编译，ForceModel.propagate 检测到任一 force 返回 None 时抛能力错误（显式报错，不静默回退到 Python eom 路径）。子类按需覆盖。元组协议见 parse_force_tuple （Rust lib.rs）：

• GravityField: ("gravity", c_flat, s_flat, mu, radius, degree, order, input_frame, propagation_frame, body, propagation_origin, tide_mode, k_love_flat, k_plus_flat_or_none)

• ThirdBody: ("third_body", naif_id_str, mu)

• Indirect: ("indirect", naif_id_str, mu)

• SRP: ("srp", area, mass, cr, shadow_bodies_list)

参数:

system (System) -- 当前动力学系统（用于查 origin / frame 等运行时参数）。

返回:

力元组，或 None。

返回类型:

tuple | None

e2m2e.algorithm.forces.physical_model.require_inertial_frame(system, t)

校验参考系为惯性系，返回 (coordinate_system, spice, origin_body)。

供在传播惯性系（ICRF，轴旋转矩阵为单位阵）中直接计算的力模型调用。 非惯性系（如 ITRFApproxAxes）抛 NotImplementedError。

参数:

• system (Any)

• t (float)

返回类型:

tuple[Any, Any, str]

球谐重力场力模型。

class e2m2e.algorithm.forces.gravity_field.GravityField(body, degree=2, order=None, gravity_file=None, input_frame=None, tide_mode='none', tide_convention='tide_free', epoch=None, polar_motion_provider=None)

基类：PhysicalModel

球谐重力场模型。

在指定的固连坐标系（默认 ITRF93）中展开球谐级数，计算引力加速度。 加速度计算全部由 Rust 编译路径承载（("gravity", ...) 力元组， crates/e2m2e-forces/src/forces/gravity_field.rs，含潮汐），Python 侧 不保留参考实现（issue #378）。

参数:

• body (str)

• degree (int)

• order (int | None)

• gravity_file (str | Path | None)

• input_frame (str | None)

• tide_mode (str)

• tide_convention (str)

• epoch (float | None)

• polar_motion_provider (Callable[[float], tuple[float, float]] | None)

__init__(body, degree=2, order=None, gravity_file=None, input_frame=None, tide_mode='none', tide_convention='tide_free', epoch=None, polar_motion_provider=None)

初始化 GravityField。

参数:

• body (str) -- 中心天体名称,如 'EARTH'、'MOON'。

• degree (int) -- 最大 degree,默认 2。

• order (int | None) -- 最大 order,默认等于 degree。

• gravity_file (str | Path | None) -- 自定义重力场文件路径（.gfc 或 .cof）；None 时按 body 取包内默认文件（地球 EGM96-to10，月球 GRGM900C）。

• input_frame (str | None) -- 球谐展开坐标系的 SPICE frame 名。None 时按 body 推导：地球 ITRF93、月球 MOON_PA；其它天体需显式提供。

• tide_mode (str) -- 潮汐档位，对齐 GMAT ETide 三档： "none" （无潮汐）、"solid" （固体潮 Step1+Step2）、 "solid_and_pole" （固体潮 + 极潮）。

• tide_convention (str) -- 系数约定，"tide_free" 或 "zero_tide" 。 zero_tide 模式减去永久潮汐（系数已含永久分量）。

• epoch (float | None) -- dot 项（系数长期变化率）外推的参考历元（SPICE et 秒）。 与 .gfc 的 dot 行配合；None 表示不外推。

• polar_motion_provider (Callable[[float], tuple[float, float]] | None) -- 极潮 xp/yp 提供者，签名 (et) -> (xp, yp) （arcsec）。solid_and_pole 档必需；由调用方从 gmat_eop 注入。

返回类型:

None

property body: str

中心天体名称。

property degree: int

最大 degree。

property order: int

最大 order。

property input_frame: str

球谐展开坐标系的 SPICE frame 名。

property gravity_file: str | Path | None

用户传入的自定义 .gfc 路径；None 表示用包内默认 EGM96。

property gravitational_parameter: float

引力参数 GM。

property reference_radius: float

参考半径 R_e。

property coefficients: dict[str, ndarray[tuple[Any, ...], dtype[floating]]]

正规化系数副本。

property tide_mode: str

潮汐档位。

to_rust_spec(system)

序列化为 Rust propagate_compiled 的 ("gravity", ...) 元组。

SolidAndPole 档暂不支持（需外部 xp/yp provider），返回 None 让 ForceModel 回退 Python 路径。

参数:

system (Any)

返回类型:

tuple | None

property tide_convention: str

系数约定。

大气阻力力模型。

class e2m2e.algorithm.forces.drag.DragModel(atmosphere, area, mass, body='EARTH', cd=2.2)

基类：PhysicalModel

大气阻力力模型。

在 ITRF（地固系）中计算大气密度与相对速度，求得阻力加速度后转换回 参考系。大气在 ITRF 中静止，因此相对速度等于航天器 ITRF 速度。

加速度计算全部由 Rust 编译路径承载（("drag", ...) 力元组， crates/e2m2e-forces/src/forces/drag.rs），Python 侧不保留参考实现 （issue #378）。to_rust_spec 需 system 提供 SPICE（ITRF93 pxform 帧旋转）；不满足时返回 None，ForceModel.propagate 据此显式报 能力错误（不静默回退）。

参数:

• atmosphere (ExponentialAtmosphere) -- 大气密度模型（依赖注入）。

• body (str) -- 中心天体名称，默认 'EARTH'。

• cd (float) -- 阻力系数，默认 2.2。

• area (float) -- 航天器迎风截面积，单位 m²。

• mass (float) -- 航天器质量，单位 kg。

property atmosphere: ExponentialAtmosphere

大气密度模型。

property body: str

中心天体名称。

property cd: float

阻力系数 Cd。

property area: float

迎风截面积，单位 m²。

property mass: float

航天器质量，单位 kg。

property ballistic_coefficient: float

弹道系数 Cd·A/m，单位 m²/kg。

to_rust_spec(system)

序列化为 Rust ("drag", area, mass, cd, propagation_frame, f107, ap) 元组。

f107/ap 从注入的大气模型取出，确保 Rust 路径与配置用同一组太阳活动 参数（issue #315 的 drag 静默分歧先例，Rust 与配置同源）。

需要 system 提供 SPICE 以做 ITRF93 pxform 帧旋转。若 system 未暴露 spice 属性、或中心天体非 EARTH，返回 None——由 ForceModel.propagate 显式报能力错误，不静默回退 Python 路径。

参数:

system (Any)

返回类型:

tuple | None

推力与机动模型。

提供两种推力/机动表示：

• ImpulsiveBurn：瞬时 Δv 机动事件，由 ForceModel.propagate_maneuvers 在指定 epoch 处中断传播并施加速度增量。

• FiniteBurn：连续推力加速度力模型，继承 PhysicalModel， 在传播过程中实时参与加速度计算。

FiniteBurn 合并了 GMAT R2026a 的 FiniteBurn （配置）与 FiniteThrust （力模型）两层，未引入 Thruster 硬件层。 推力大小与方向解耦：thrust_profile(t) 返回标量推力（N）， direction 给出方向向量（固定向量或随状态更新的可调用）， 内部归一化为单位向量。质量为常量（不支持推进剂消耗）。

VariableMassFiniteBurn 是其可变质量对应物：质量作为状态量 state[6] 随推力消耗（ṁ = −T/(Isp·g₀)），是低推力最优控制与 月面动力下降的受控动力学基座。详见 docs/plans/lowthrust-foundation-prd.md。

direction_frame 支持 "VNB"、"LVLH" 与 None：

• None：direction 直接在传播（惯性）坐标系内解释。

• "VNB"：direction 在 VNB 坐标系下解释，其中 \(V = v/\\|v\\|\) (速度方向)， \(N = (r \\times v)/\\|r \\times v\\|\) (角动量方向)， \(B = V \\times N\) (副法向)。

• "LVLH"：direction 在 LVLH 坐标系下解释，其中 \(R = r/\\|r\\|\) (径向)， \(N = (r \\times v)/\\|r \\times v\\|\) (法向)， \(T = N \\times R\) (沿迹方向)。

class e2m2e.algorithm.forces.thrust.ImpulsiveBurn(epoch, delta_v)

基类：object

瞬时 Δv 机动事件。

delta_v 在传播（惯性）坐标系内解释，由 e2m2e.algorithm.forces.force_model.ForceModel.propagate_maneuvers() 在 epoch 处施加 state[3:6] += delta_v。

VNB/LVLH burn 坐标系暂不支持（届时加 frame 字段，转换走 CoordinateSystem.transform_vector()，对应 GMAT Burn::ConvertDeltaVToInertial 的 coincident=true 纯旋转）。

参数:

• epoch (float) -- 施加时刻，SPICE et 秒，与 ForceModel.propagate 的 t_span 一致。

• delta_v (ndarray[tuple[Any, ...], dtype[floating]]) -- 速度增量，参考系，形状 (3,)。

epoch: float

delta_v: ndarray[tuple[Any, ...], dtype[floating]]

class e2m2e.algorithm.forces.thrust.FiniteBurn(thrust_profile, direction, mass, direction_frame=None)

基类：PhysicalModel

恒质量连续推力加速度力模型。

6D 状态传播由 Rust 编译路径执行。配置 DSL 构造的常量或 pulse 推力曲线和 固定方向可下沉；任意 Python callable 无法进入 Rust RK 内循环，在传播入口会 显式报能力错误。需要推进剂消耗时使用 VariableMassFiniteBurn（变质量，7D 状态）。

direction 给出方向向量（固定向量或随状态更新的可调用）， 内部归一化为单位向量。质量为常量（不支持推进剂消耗）。

direction_frame 支持 "VNB"、"LVLH" 与 None：

• None：direction 直接在传播（惯性）坐标系内解释。

• "VNB"：三个分量依次对应速度单位向量、角动量单位向量和副法向量。

• "LVLH"：三个分量依次对应径向单位向量、沿迹单位向量和轨道面法向量； 沿迹单位向量由法向量叉径向量得到。

参数:

• thrust_profile (Callable[[float], float]) -- t -> thrust （N，标量；0 表示关机）。

• direction (npt.ArrayLike | Callable[[float, npt.NDArray[np.floating]], npt.ArrayLike]) -- 固定方向向量 (3,)，或 (t, state) -> (3,) 可调用。

• mass (float) -- 航天器质量（kg，常量）。

• direction_frame (str | None) -- 方向解释坐标系，"VNB" / "LVLH" / None。

property thrust_profile: Callable[[float], float]

推力大小随时间变化的可调用（N）。

property direction: ArrayLike | Callable[[float, ndarray[tuple[Any, ...], dtype[floating]]], ArrayLike]

推力方向：固定向量或 (t, state) -> (3,) 可调用。

property direction_frame: str | None

方向解释坐标系：'VNB'、'LVLH' 或 None。

property mass: float

航天器质量（kg，常量）。

to_rust_spec(system)

序列化为恒质量 6D 编译传播接受的推力规格。

只有配置 DSL 构造的 constant/pulse 推力 profile 和固定方向可以下沉 到 Rust；任意 Python callable 无法在 Rust RK 内安全求值，返回 None。 返回规格为 ("low_thrust", mass, thrust, t_start, t_end, direction, direction_frame)，其中 constant profile 的起止时间为 None。

参数:

system (object)

返回类型:

tuple | None

class e2m2e.algorithm.forces.thrust.VariableMassFiniteBurn(thrust, isp, initial_mass, direction, direction_frame=None)

基类：PhysicalModel

连续推力加速度力模型（质量随推力消耗）。

与 FiniteBurn 的唯一区别：质量不是常量，而是状态量 state[6]。低推力转移与月面动力下降等最优控制问题中，质量演化 是燃耗最优的基本变量（ṁ = −T/(Isp·g₀)），必须纳入状态向量。

配套的 7D 传播在 propagate 中走 Rust 快速路径 propagate_compiled_lowthrust：状态 [x, y, z, vx, vy, vz, m]，受控动力学在 Rust 侧（augmented_state 的 augmented_eom_7d）。详见 docs/plans/lowthrust-foundation-prd.md。

推力大小与方向解耦，语义同 FiniteBurn。direction 支持 固定向量或 (t, state) -> (3,) 可调用；state 为 7D 时可调用方向 可读取 state[6] 中的质量。direction_frame 支持 "VNB" / "LVLH" / None，帧解析与 FiniteBurn 一致。

参数:

• thrust (float) -- 推力幅值（N，常量）。

• isp (float) -- 比冲（s）。

• initial_mass (float) -- 初始质量（kg），用于初始化状态第 7 维与校验。

• direction (npt.ArrayLike | Callable[[float, npt.NDArray[np.floating]], npt.ArrayLike]) -- 固定方向向量 (3,)，或 (t, state) -> (3,) 可调用。

• direction_frame (str | None) -- 方向解释坐标系，"VNB" / "LVLH" / None。

property thrust: float

推力幅值（N，常量）。

property isp: float

比冲（s）。

property initial_mass: float

初始质量（kg），用于初始化状态第 7 维。

property direction: ArrayLike | Callable[[float, ndarray[tuple[Any, ...], dtype[floating]]], ArrayLike]

推力方向：固定向量或 (t, state) -> (3,) 可调用。

property direction_frame: str | None

方向解释坐标系：'VNB'、'LVLH' 或 None。

to_rust_spec(system)

序列化为低推力 7D 传播路径接受的推力规格。

仅当 direction 为固定向量时返回元组（可调用方向需 Python 求值， 无法下沉到 Rust）；常量推力映射成满油门（throttle = 1.0）， t_max = thrust。返回元组会被 ForceModel 的低推力分支拆出， 交给 propagate_compiled_lowthrust，不经过 6D 的 CompiledForce 路径。

参数:

system (object)

返回类型:

tuple | None

class e2m2e.algorithm.forces.thrust.BurnApplication(index, epoch, delta_v, velocity_before, velocity_after)

基类：object

单次脉冲机动在 propagate_maneuvers 输出中的记录。

参数:

• index (int) -- post-burn 行在输出 states 中的行号。

• epoch (float) -- 施加时刻，SPICE et 秒。

• delta_v (ndarray[tuple[Any, ...], dtype[floating]]) -- 施加的 Δv，传播坐标系，(3,)。

• velocity_before (ndarray[tuple[Any, ...], dtype[floating]]) -- 施加前速度，(3,)。

• velocity_after (ndarray[tuple[Any, ...], dtype[floating]]) -- 施加后速度，(3,)。

index: int

epoch: float

delta_v: ndarray[tuple[Any, ...], dtype[floating]]

velocity_before: ndarray[tuple[Any, ...], dtype[floating]]

velocity_after: ndarray[tuple[Any, ...], dtype[floating]]

第三体引力间接项。

class e2m2e.algorithm.forces.indirect_term.IndirectTerm(body, mu=None)

基类：PhysicalModel

第三体引力的间接项（geocentric 加速系所需）。

在以某天体（如地球）为原点的非惯性系下传播时，运动方程需对每个摄动 天体 \(i\) 补一项 -μ_i · r_i / |r_i|³ （间接项），扣除摄动天体 对原点的引力（见 EphemerisDynamics 的 N 体闭式公式）。

ThirdBodyGravity 内部已自带间接项，但 GravityField 只算球谐 直接引力（含中心项 degree=0），不带间接项。所以用 GravityField 模拟月球（中心+非球形）时，必须单独补月球间接项——既不能用 ThirdBodyGravity("MOON") （会与 GravityField 的 degree=0 中心项 重复算月球点质量），也不能省略（地心系下物理不正确）。

加速度：-μ_body · r_body / |r_body|³，其中 r_body 为摄动天体相对 system.origin 的位置（由 system.get_body_position 自动以 origin 为观察者计算）。与 ThirdBodyGravity 的间接项逐字一致。

参数:

• body (str) -- 摄动天体名称（如 'MOON'）。

• mu (float | None) -- 引力参数（km³/s²）。为 None 时， 在 to_rust_spec 中从 system.gravitational_parameter(body) 获取。

property body: str

摄动天体名称。

property mu: float | None

显式设置的引力参数；None 表示从 system 获取。

to_rust_spec(system)

序列化为 ("indirect", naif_id_str, mu)。

返回类型:

tuple | None

相对论修正力模型。

class e2m2e.algorithm.forces.relativistic_correction.RelativisticCorrection(central_body, *, primary_body='SUN', enable_schwarzschild=True, enable_lense_thirring=True, enable_de_sitter=True, angular_momentum_vector=None, body_radius=None, c=299792.458, gamma=1.0)

基类：PhysicalModel

相对论修正力模型。

实现 Schwarzschild、Lense-Thirring 与 de Sitter（geodesic）三项相对论 加速度修正，公式与 GMAT R2026a RelativisticCorrection 对齐。

加速度计算全部由 Rust 编译路径承载（("relativistic", ...) 力元组， crates/e2m2e-forces/src/forces/relativistic.rs），Python 侧不保留参考 实现（issue #378）。

参数:

• central_body (str)

• primary_body (str | None)

• enable_schwarzschild (bool)

• enable_lense_thirring (bool)

• enable_de_sitter (bool)

• angular_momentum_vector (npt.ArrayLike | None)

• body_radius (float | None)

• c (float)

• gamma (float)

property central_body: str

中心天体名称（大写）。

property primary_body: str | None

de Sitter 项主天体名称（大写），可能为 None。

property enable_schwarzschild: bool

Schwarzschild 项开关。

property enable_lense_thirring: bool

Lense-Thirring 项开关。

property enable_de_sitter: bool

de Sitter 项开关。

property angular_momentum_vector: ndarray[tuple[Any, ...], dtype[floating]] | None

Lense-Thirring 角动量矢量（覆盖值），单位 km²/s。

注意：这里的 J 与 GMAT 约定一致，是 (2/5) * R² * spin_rate 形式的归一化量，不是 SI 物理角动量（kg·m²/s）。

property body_radius: float | None

中心天体赤道半径（覆盖值），单位 km。

property c: float

光速，单位 km/s。

property gamma: float

后牛顿参数 gamma。

to_rust_spec(system)

序列化为 ("relativistic", ...) 元组。

• LT 项需要 sxform + body-fixed frame；本仓库已实测 NRHO 上 LT 量级 < 1m （#343 排查），但完整移植已实现（含 sxform via cspice-sys FFI）。

• 如果 LT 启用但 angular_momentum_vector 未传，Rust 侧会每步 sxform 自动算（与 Python 一致）；如需避免 sxform 开销，可在 Python 侧 预先算好 J 向量并传入 angular_momentum_vector。

返回类型:

tuple | None

太阳辐射压与阴影模型

e2m2e 提供基于 cannonball 模型的太阳辐射压（SRP）力模型，以及圆锥阴影模型用于计算地影/月影对光照的遮挡效应。两者均通过 PhysicalModel 接口与力模型组合集成，支持配置驱动的序列化与反序列化。

太阳辐射压模型

SolarRadiationPressure 实现 Montenbruck & Gill 的 cannonball SRP 模型：

\[\mathbf{a} = \text{flux} \cdot P_{1\text{AU}} \left(\frac{1\ \text{AU}}{r}\right)^2 \frac{C_R \, A}{m} \, \hat{\mathbf{u}}\]

其中 P_1AU = 4.56e-6 N/m² 为 1 AU 处太阳光压常数，r 为航天器到太阳的距离， C_R 为辐射反射系数（1 = 全吸收，2 = 全反射），A 为迎风截面积（m²）， m 为质量（kg）。flux ∈ [0, 1] 由阴影模型给出，全光照为 1，本影为 0。

参数说明：

参数	含义	默认值
area	航天器迎风截面积（m²）	必填
mass	航天器质量（kg）	必填
cr	辐射反射系数 C_R	1.5
shadow	阴影模型实例（注入）	None （全光照）

阴影模型

ConicalShadowModel 定义 flux_factor(t, state, system) -> float 接口，是当前唯一的阴影模型实现。

圆锥阴影模型

实现 GMAT ShadowState 的圆锥阴影算法（Montenbruck & Gill §3.4.2）， 从航天器看太阳与遮挡体的视角径 (a, b) 与角距 c，分四分支判定：

• 全光照：遮挡体与太阳圆盘不相交

• 本影：遮挡体完全遮住太阳圆盘（flux = 0）

• 半影：部分重叠，用 M&G eq. 3.92-3.94 精确圆面重叠面积计算

• 环形食：遮挡体小于太阳，中心对齐但边缘透光

多遮挡体（如地球 + 月球）的光照份额合成遵循 GMAT GMT-6543 规范： 任一遮挡体本影 → 0；两体部分阴影且不重叠 → 包容排斥；重叠 → 保守取最小值。

参数说明：

参数	含义	默认值
bodies	遮挡体名称列表（大写）	("EARTH",)
radii	天体半径覆盖字典（km）	None （使用内置默认值）

内置默认半径：

天体	半径（km）
EARTH	6378.1363
MOON	1737.4
SUN	695700.0

配置与序列化

SRP 与阴影模型支持通过 force_config 进行配置驱动的序列化。

配置字典格式：

{
    "type": "SolarRadiationPressure",
    "params": {
        "area": 2.0,
        "mass": 1000.0,
        "cr": 1.5,
        "shadow": {
            "type": "ConicalShadowModel",
            "params": {
                "bodies": ["EARTH", "MOON"],
                "radii": None
            }
        }
    }
}

SRP 工作流示例

以下示例展示 SRP + 地影/月影 + EphemerisDynamics 的传播流程：

import numpy as np

from e2m2e.algorithm.coordinate import (
    CelestialBodyOrigin,
    CoordinateSystem,
    ICRSAxes,
)
from e2m2e.algorithm.dynamics import EphemerisSystem
from e2m2e.data.kernels.manager import SPICEManager
from e2m2e.algorithm.forces import (
    SolarRadiationPressure,
    ConicalShadowModel,
    ForceModel,
)

# 1. 加载 SPICE 内核
mgr = SPICEManager()
mgr.load_kernel("path/to/de440.bsp")

# 2. 构建星历系统（含地球、月球、太阳；frame 默认 J2000）
system = EphemerisSystem(
    bodies=["EARTH", "MOON", "SUN"],
    spice=mgr,
    origin="EARTH",
)

# 3. 设置坐标系（ICRF 惯性系，力模型要求惯性系）
axes = ICRSAxes()
origin = CelestialBodyOrigin(body="EARTH", spice=mgr)
system.coordinate_system = CoordinateSystem(axes=axes, origin=origin)

# 4. 创建阴影模型（地影 + 月影）
shadow = ConicalShadowModel(bodies=["EARTH", "MOON"])

# 5. 创建 SRP 力模型
srp = SolarRadiationPressure(
    area=2.0,      # m²
    mass=1000.0,   # kg
    cr=1.5,
    shadow=shadow,
)

# 6. 组装 ForceModel
fm = ForceModel(system)
fm.add_force(srp, name="SRP")

# 7. 初始状态：LEO 近似圆轨道
r0 = 6678.0  # km（约 300 km 高度）
v0 = np.sqrt(398600.435507 / r0)  # km/s
state0 = np.array([r0, 0.0, 0.0, 0.0, v0, 0.0])

# 8. 传播（1 个轨道周期，约 90 分钟）
et0 = mgr.utc_to_et("2024-06-21T00:00:00")
period = 2 * np.pi * r0 / v0  # ~5400 s
result = fm.propagate(state0, (et0, et0 + period))

print(f"传播点数: {len(result['time'])}")
print(f"末状态: {result['states'][-1]}")

# 9. 序列化配置到 JSON
from e2m2e.algorithm.forces import dump_force_config
dump_force_config(fm, "srp_config.json")

# 10. 从 JSON 恢复
from e2m2e.algorithm.forces import load_force_config
fm2 = load_force_config("srp_config.json", system)

# 验证 round-trip
assert fm2.to_config() == fm.to_config()

# 11. 启用/禁用 SRP 对比
fm.disable("SRP")
result_no_srp = fm.propagate(state0, (et0, et0 + period))

fm.enable("SRP")
result_with_srp = fm.propagate(state0, (et0, et0 + period))

# 对比末位置差异
diff = np.linalg.norm(result_with_srp["states"][-1, :3]
                     - result_no_srp["states"][-1, :3])
print(f"SRP 引起的 1 周期位置差异: {diff:.3f} km")

纯函数测试路径

SRP 和阴影模型均提供纯函数接口，可在无 SPICE 环境下直接测试：

import numpy as np
from e2m2e.algorithm.forces.srp import SolarRadiationPressure
from e2m2e.algorithm.forces.shadow import ConicalShadowModel

# SRP 纯函数测试
srp = SolarRadiationPressure(area=2.0, mass=1000.0, cr=1.5)
sun_to_sc = np.array([1.0, 0.0, 0.0]) * 149597870.691  # 1 AU
accel = srp._compute_srp_acceleration(sun_to_sc, flux_factor=1.0)
print(f"1 AU 处全光照 SRP 加速度: {accel} km/s²")

# 阴影模型纯函数测试
shadow = ConicalShadowModel()
sc_pos = np.array([7000.0, 0.0, 0.0])
body_pos = np.array([0.0, 0.0, 0.0])
sun_pos = np.array([1.5e8, 0.0, 0.0])
flux = shadow._body_flux_factor(
    sc_pos, body_pos, sun_pos,
    body_radius=6378.1363, sun_radius=695700.0
)
print(f"地影光照份额: {flux}")