Universal ODE Kernel

> Aligned with PCR Master Blueprint v1.0 — see Blueprint §2.5（dynamics_core: Universal Kernel）、§7.16（dynamics_core 严格普适）、§4.1（积分阶段归位）。 > 职责：本文定义 dynamics_core/ 的 universal kernel：物理状态原语（SpatialState / InertialState / AuxState）、运动方程（compute_state_derivative）、积分器（integrate_rk4 / integrate_euler）。这套机制不知道火箭长什么样——它只接收 Forces 与 Inertia，输出下一时刻状态。 > 核心约束：dynamics_core/ 的所有文件，grep -r "plant/\|avionics/\|fcc/\|bus/\|environment/\|simulation/" 必须零命中。 > C-Distillation note：SpatialState 是 PoD；Integrator 是纯函数；StateVector 蒸馏后是 double[16]，pack/unpack 是 memcpy。

1. Universal Kernel 的核心承诺

dynamics_core/ 是积分引擎，舞台换了，演员和导演不变（Blueprint §1）：把地球换月球、火箭换飞机，dynamics_core/ 一字不改。它的职责被严格限制为：

职责	不属于此处的
数值积分（RK4 / Euler）	❌ 力的本构计算（→ `plant/physics/`）
运动方程（牛顿 + 欧拉）	❌ 气动 / 推力 / 重力模型（→ `plant/` + `environment/`）
物理状态原语类型	❌ 火箭硬件 / 设备（→ `plant/model/` + `avionics/`）
拓扑代数（`StageOp`）	❌ 分离事件来源（→ `simulation::interpret_event`）
Forces Monoid 类型	❌ Forces 的具体填充（→ `plant/physics/*` 各力源）

> 关键测试：dynamics_core/ 任何头文件不得 #include plant/、avionics/、fcc/、bus/、environment/、simulation/。

2. 目录与依赖

src/dynamics_core/
├── state/
│   ├── SpatialState.h          ← 位姿原语（pos_lic / vel_lic / quat / omega_body）
│   ├── InertialState.h         ← 质量+转动惯量
│   ├── AuxState.h              ← IMU 真值瞬态（spec_force_b + omega_b）
│   └── KinematicsView.h        ← 强类型 view 模板
├── ode/
│   ├── EquationsOfMotion.h     ← 纯函数导数计算
│   ├── EquationsOfMotion.cpp
│   ├── Integrator.h            ← integrate_rk4 / integrate_euler
│   └── Integrator.cpp
├── pipeline/
│   ├── Forces.h                ← Monoid 类型（详见 Forces_Monoid.md）
│   └── RWSAlias.h              ← BodyRWS / WorldRWS 模板别名
└── algebra/
    ├── StageOp.h               ← WorldStage / StageOp（详见 Topology_Algebra.md）
    └── evolve_topology.h

只允许的依赖：types/（Vec3 / Quat / Matrix / Frames） + frames/ + monad/ + contracts/。

3. 物理状态原语

3.1 SpatialState：位姿"单一事实来源"

cpp

// dynamics_core/state/SpatialState.h
namespace dynamics {

class SpatialState {
private:
    Vec3_T&lt;frame::LIC&gt;   pos_lic_;            // 发射惯性系位置
    Vec3_T&lt;frame::LIC&gt;   vel_lic_;            // 发射惯性系速度
    Quat                 att_lic_to_body_;    // 姿态：LIC → BODY
    Vec3_T&lt;frame::BODY&gt;  omega_body_;         // 本体系角速度

public:
    Vec3_T&lt;frame::LIC&gt;   position()         const { return pos_lic_; }
    Vec3_T&lt;frame::LIC&gt;   velocity()         const { return vel_lic_; }
    Quat                 attitude()         const { return att_lic_to_body_; }
    Vec3_T&lt;frame::BODY&gt;  angular_velocity() const { return omega_body_; }

    // 纯函数：施加导数，返回新 SpatialState（不变性）
    SpatialState apply_derivative(const SpatialDerivative& d, double dt) const;
};

}

核心约束：

状态只保存 LIC（发射惯性系）下的绝对量；所有其他视角（ECF/AERO/NUE）由 KinematicsView 在运行时从 (Reader 配置 + SpatialState) 推导
不缓存任何转换矩阵——SSOT（Single Source of Truth）
强类型 Vec3_T<frame::LIC> 防止跨系误用

3.2 InertialState：质量 + 转动惯量

cpp

// dynamics_core/state/InertialState.h
struct InertialState {
    double  mass;        // 总质量
    Matrix  inertia;     // 体系转动惯量张量（3x3，约定本体系）
    double  fuel_mass;
    double  oxidizer_mass;
};

3.3 AuxState：IMU 真值瞬态

cpp

// dynamics_core/state/AuxState.h
struct AuxInstant {
    Vec3_T&lt;frame::BODY&gt;  spec_force_b;   // 比力 = (Forces.total_force_b - mass*g_body) / mass
    Vec3_T&lt;frame::BODY&gt;  omega_b;        // 角速度（snapshot）
};

> AuxInstant 是积分计算的副产物，专门给 IMU 真值模型读（avionics/imu/ 用它产生加表/陀螺真值脉冲）。它是 dynamics_core 内部产出的"瞬态量"，不参与状态向量。

3.4 StateVector + StateMap：积分用的扁平打包

cpp

// dynamics_core/ode/EquationsOfMotion.h
using StateVector = std::vector&lt;double&gt;;  // 蒸馏后退化为 double[16]

struct StateMap {
    static constexpr size_t POS_START   = 0;   // 3 doubles
    static constexpr size_t VEL_START   = 3;
    static constexpr size_t ATT_START   = 6;   // 4 doubles (quat)
    static constexpr size_t OMEGA_START = 10;
    static constexpr size_t MASS_IDX    = 13;
    static constexpr size_t FUEL_IDX    = 14;
    static constexpr size_t OX_IDX      = 15;
    static constexpr size_t STATE_SIZE  = 16;
};

StateVector pack_state(const SpatialState& s, const InertialState& i);
void        unpack_state(const StateVector& v, SpatialState& s, InertialState& i);

> StateVector 是积分器内部的 ABI——RK4 需要 +/* 的扁平表示。pack/unpack 是 O(16) memcpy，性能可忽略。蒸馏后是栈分配 double[16]。

4. 运动方程（纯函数）

4.1 顶层接口

cpp

// dynamics_core/ode/EquationsOfMotion.h
StateVector compute_state_derivative(
    const StateVector& state,
    const Forces&      forces,        // Forces Monoid 累加结果（已在外部完成）
    const Matrix&      inertia        // 转动惯量
);

契约：

输入：当前扁平状态 + 已聚合的 Forces + 当前转动惯量
输出：每个分量的时间导数（同样扁平）
不知道力从哪来——Forces 已经在外部被 compute_thrust + compute_aero + compute_gravity 累加完毕
不知道火箭长什么样——只看 mass / inertia / spatial

4.2 分量导数（私有，纯函数）

cpp

// 速度导数：力学第二定律 + 重力（LIC 系直接累加）
Vec3 compute_velocity_derivative(const Vec3& force_b,
                                  double mass,
                                  const Quat& att_lic_to_body,
                                  const Vec3& gravity_lic) {
    Vec3 spec_force_b   = force_b / mass;
    Vec3 spec_force_lic = att_lic_to_body.conjugate().rotate(spec_force_b);
    return spec_force_lic + gravity_lic;
}

// 位置导数：直接等于速度
Vec3 compute_position_derivative(const Vec3& velocity_lic) { return velocity_lic; }

// 角速度导数：欧拉方程 I·ω̇ = M − ω × (I·ω)
Vec3 compute_angular_velocity_derivative(const Matrix& I,
                                          const Vec3& omega_b,
                                          const Vec3& moment_b) {
    Vec3 gyro_torque = omega_b.cross(I * omega_b);
    return I.inverse() * (moment_b - gyro_torque);
}

// 四元数导数：q̇ = ½·q·(0, ω)
Quat compute_attitude_derivative(const Quat& q, const Vec3& omega_b) {
    Quat w_quat(0, omega_b.x(), omega_b.y(), omega_b.z());
    Quat dq = q * w_quat;
    dq.scale(0.5);
    return dq;
}

重力的特殊通道：注意 gravity_lic 是 Forces 里独立的字段（详见 Forces_Monoid.md §3）——重力直接以 LIC 加速度形式累加，不进 total_force_b。这是因为 gravity_lic 在 LIC 系是天然的标量场表示，强行投到 BODY 再投回会引入数值误差。

5. 积分器：可插拔策略

5.1 函数式签名

cpp

// dynamics_core/ode/Integrator.h
namespace dynamics::ode {

using IntegratorFunc = std::function&lt;
    StateVector(const StateVector&  /*y*/,
                const Forces&        /*forces*/,
                const Matrix&        /*inertia*/,
                double               /*dt*/)
&gt;;

// 一阶欧拉
StateVector integrate_euler(const StateVector& y,
                             const Forces&,
                             const Matrix& inertia,
                             double dt);

// 四阶 Runge-Kutta（默认）
StateVector integrate_rk4(const StateVector& y,
                           const Forces&,
                           const Matrix& inertia,
                           double dt);

}

5.2 RK4 实现要点

cpp

StateVector integrate_rk4(const StateVector& y,
                           const Forces& forces,
                           const Matrix& inertia,
                           double dt) {
    auto k1 = compute_state_derivative(y, forces, inertia);
    auto k2 = compute_state_derivative(y + (dt/2) * k1, forces, inertia);
    auto k3 = compute_state_derivative(y + (dt/2) * k2, forces, inertia);
    auto k4 = compute_state_derivative(y + dt * k3, forces, inertia);
    return y + (dt / 6) * (k1 + 2*k2 + 2*k3 + k4);
}

简化前提：Forces 与 inertia 在 [t, t+dt] 内视作常量。这是一个工程取舍——严格 RK4 需要每个子步重算力，但火箭场景下 dt=10ms 内力变化平缓，常量近似的误差远小于本构模型本身的精度。

> 如果某个 mission 需要 sub-step force 重算，未来可加 IntegratorFuncWithForcesFn——签名变为 (y, force_fn, inertia, dt)，其中 force_fn 是 (y) → Forces 闭包。这是非默认路径。

5.3 单步积分的 RWS 包装（位于 simulation/，非 dynamics_core/）

> ⚠️ 重要归属：把 integrate_* 嵌入 BodyRWS pipeline 的代码不在 dynamics_core/ —— 它在 simulation/pipeline/BodyTick.cpp 里，因为 BodyRWS pipeline 跨域（Blueprint §2.6.4）。dynamics_core/ 只暴露纯函数 integrate_rk4(y, forces, inertia, dt)。

cpp

// simulation/pipeline/BodyTick.cpp（示意，归属 simulation/，本文档不展开）
BodyRWS&lt;std::monostate&gt; integrate_step(IntegratorFunc fn, double dt) {
    return body_get() &gt;&gt;= [fn, dt](RocketBody body) {
        // 拿到 Forces（来自上游 plant/physics 累加）
        // pack → integrate → unpack → 更新 AuxState
        // ...
        return body_put(updated);
    };
}

详见 06_Simulation/Body_World_Tick.md。

6. KinematicsView：从 LIC SSOT 投影到其他坐标系

cpp

// dynamics_core/state/KinematicsView.h
template &lt;typename TargetFrame&gt;
struct KinematicsView {
    Vec3_T&lt;TargetFrame&gt; position;
    Vec3_T&lt;TargetFrame&gt; velocity;
    Quat                attitude_to_body;   // TargetFrame → BODY
};

KinematicsView 是纯结构——dynamics_core/ 只定义类型，不定义如何从 SpatialState 提取（提取需要时间、发射点 ECF 坐标等，这些是 BodyEnv / WorldEnv 的责任，属于 simulation/）。

cpp

// simulation/env/TransformContext.h（不在 dynamics_core 内）
struct EcfTransformContext {
    Vec3 launch_ecf;
    Quat lic_to_ecf;
    Time t;
    KinematicsView&lt;frame::ECF&gt; extract_view(const SpatialState& s) const;
};

> dynamics_core/ 提供模板与原语；simulation/ 提供具体的 context provider。两者通过 KinematicsView<F> 类型互通。

7. 状态推进（apply_derivative）

cpp

// dynamics_core/state/SpatialState.cpp
SpatialState SpatialState::apply_derivative(const SpatialDerivative& d, double dt) const {
    return SpatialState(
        pos_lic_         + d.d_position * dt,
        vel_lic_         + d.d_velocity * dt,
        (att_lic_to_body_ + d.d_attitude * dt).normalized(),   // 四元数归一化
        omega_body_      + d.d_omega * dt
    );
}

SpatialDerivative：

cpp

struct SpatialDerivative {
    Vec3_T&lt;frame::LIC&gt;   d_position;
    Vec3_T&lt;frame::LIC&gt;   d_velocity;
    Quat                 d_attitude;      // 注意：quat 加法是分量加，不是旋转复合
    Vec3_T&lt;frame::BODY&gt;  d_omega;
};

四元数加法+归一化是 Euler 积分的标准技巧；RK4 内部多次 derivative 累加后最终再 normalize 一次即可。

8. 与上下游的契约

8.1 上游：Forces 累加器（来自 plant/physics/）

text

plant::physics::compute_thrust_contribution(...)  → Forces (推力部分)
plant::physics::compute_aero_contribution(...)    → Forces (气动部分)
plant::physics::compute_gravity_contribution(...) → Forces (gravity_lic 通道)
                            │
                            ▼
                     Forces Monoid + 累加（详见 Forces_Monoid.md §3）
                            │
                            ▼
            compute_state_derivative(state, forces, inertia)
                            │
                            ▼
                       integrate_rk4 (此 kernel)
                            │
                            ▼
                       unpack_state → SpatialState + InertialState

8.2 下游：AuxInstant（给 IMU 真值用）

cpp

plant::AuxInstant compute_aux_instant(const StateVector& state, const Forces& forces);

> ⚠️ 实际代码现位置：src/dynamics/ode/EquationsOfMotion.h 暴露 plant::AuxInstant 返回类型——这是 v0 遗留违规（dynamics_core 不应引用 plant）。v1 迁移：把 AuxInstant 移到 dynamics_core/state/AuxState.h（已在 §3.3 列出位置），plant/ 引用 dynamics 类型而不是反过来。

8.3 不耦合：dynamics_core 不知道的事情

dynamics_core 看不见	提供者	通道
火箭长啥样	`plant/model/`	不通过任何接口；plant/ 自己用 `Forces` 给 dynamics 喂值
气动表 / 推力曲线	`plant/assets/`	同上
Atmosphere / Gravity Field	`environment/`	同上
FCC 决策	`fcc/`	完全不感知
Bus 协议	`bus/`	完全不感知
多体协调	`simulation/pipeline/WorldTick`	dynamics_core 只看单体导数

9. 反模式

反模式	为什么不行
`dynamics_core/ode/Integrator.cpp` 内 `#include "plant/..."`	违反 universal 承诺（Blueprint §2.5 / §7.16）；grep 测试会红
把 `Atmosphere` / `EngineCurve` 塞进 `Forces`	Forces 只是数值容器；本构计算在 plant/
Integrator 内部硬编码 `if (stage == COMPOSITE_FLIGHT)`	dynamics_core 不该知道 WorldStage 决定何种导数；用 `Forces` 表达即可
`SpatialState` 缓存 ECF / NUE 转换矩阵	破坏 SSOT；用 `KinematicsView` + Context 在运行时提取
在 `dynamics_core/state/` 直接保存 `RocketBody`	RocketBody 是跨域复合体（Blueprint §2.6.1）；dynamics_core 只持有"物理状态原语"
RK4 子步内重算 Forces 但 Forces 又依赖 plant/	已是 simulation/ 范畴；dynamics_core 接受 `Forces const&` 即可
把 `pack_state/unpack_state` 暴露给 plant/ 直接用	它是积分器内部 ABI；plant/ 应通过 `SpatialState` accessor 而不是扁平向量

10. C-Distillation 视角

类型	C++ 阶段	C 蒸馏阶段
`StateVector`	`std::vector<double>`	`double[16]`（栈分配）
`Quat` / `Vec3` / `Matrix`	模板/类	`struct { double x,y,z; }` / `double m[9]`
`IntegratorFunc`	`std::function`	函数指针 `StateVector(*)(...)`
`SpatialState`	class with getters	PoD struct
`KinematicsView<F>`	模板	per-frame typedef PoD
`apply_derivative`	method	static function

整个 kernel 蒸馏后是一组无依赖的 C 自由函数 + PoD struct，可单独编译进固件 ROM。

11. 测试策略

测试目标	级别	描述
`dynamics_core/` 零业务依赖	structural	`grep -r "plant/\|avionics/\|fcc/\|bus/\|environment/\|simulation/" src/dynamics_core/` 必须空
积分器收敛	bench	自由落体 100s，RK4 解 vs 解析解误差 < 1e-6
守恒律	bench	无外力（`Forces::zero()`）下 1000s 积分，总能量漂移 < 1e-9
四元数归一化稳定	bench	持续 1Hz 角速度自旋 3600s，‖q‖ 偏离 1 < 1e-12
pack/unpack 同构	unit	`unpack(pack(s, i)) == (s, i)`
强类型防呆	static	`Vec3_T<LIC> + Vec3_T<BODY>` 编译失败

12. Cross References

Forces Monoid 详细 → Forces_Monoid.md
Topology 代数与 WorldStage → Topology_Algebra.md
BodyRWS pipeline 跨域实例化 → 06_Simulation/Body_World_Tick.md
坐标系契约 → 01_Foundation/Coordinate_Frames.md + Contracts_PerBodyFrames.md
强类型 Vec3_T → 01_Foundation/Coordinate_Frames.md §3
现有代码 → src/dynamics/{ode/, state/SpatialState.h}（待 Wave 0 迁移到 src/dynamics_core/）
Blueprint §2.5（dynamics_core 定义）、§7.16（严格普适）、§4.1（推力数据流积分阶段）

Universal ODE Kernel ​

1. Universal Kernel 的核心承诺 ​

2. 目录与依赖 ​

3. 物理状态原语 ​

3.1 SpatialState：位姿"单一事实来源" ​

3.2 InertialState：质量 + 转动惯量 ​

3.3 AuxState：IMU 真值瞬态 ​

3.4 StateVector + StateMap：积分用的扁平打包 ​

4. 运动方程（纯函数） ​

4.1 顶层接口 ​

4.2 分量导数（私有，纯函数） ​

5. 积分器：可插拔策略 ​

5.1 函数式签名 ​

5.2 RK4 实现要点 ​

5.3 单步积分的 RWS 包装（位于 simulation/，非 dynamics_core/） ​

6. KinematicsView：从 LIC SSOT 投影到其他坐标系 ​

7. 状态推进（apply_derivative） ​

8. 与上下游的契约 ​

8.1 上游：Forces 累加器（来自 plant/physics/） ​

8.2 下游：AuxInstant（给 IMU 真值用） ​

8.3 不耦合：dynamics_core 不知道的事情 ​

9. 反模式 ​

10. C-Distillation 视角 ​

11. 测试策略 ​

12. Cross References ​