# Numerical-Relativity Spacetime Movie Renderer ## 1. Repo 目标 本项目的目标是开发一个**离线、尽可能物理正确的数值相对论时空视频渲染器**，最终用于渲染双黑洞（BBH）并合等动态数值时空的 4K 视频。 核心目标不是实时渲染，也不是构造视觉上“像双黑洞”的近似度规，而是： - 直接消费数值相对论演化输出的真实 4D 时空； - 对每一帧的相机光线做完整的 time-dependent backward null ray tracing； - 支持真实巡天恒星数据作为无穷远天球背景； - 正确处理恒星点源的多像、放大、频移和最终 PSF； - 允许相机沿一般 4D worldline 运动，并使用预先生成的 tetrad/标架轨迹； - 让平直时空、解析时空、数值时空在同一渲染框架中作为可替换 backend； - 最终能够“看到”每次 NR 代码实际跑出来的时空，而不是只看 waveform 或标量诊断。 第一阶段不考虑物质辐射、吸积盘、流体、等离子体等局域发射源。每条 ray 的正常终点 为：逃逸到某个无穷远天球，或达到红移暗阈值。预算耗尽与数据/积分失败是单独的 `UNRESOLVED` / `INCOMPLETE` 类别，不与物理暗终态混同（见 §18）。 --- ## 2. 明确的非目标 当前 repo 不以以下内容为主要目标： - 实时渲染； - ShaderToy / 游戏式视觉近似； - frozen-snapshot 作为最终物理方案； - 直接渲染吸积盘、喷流或其他 matter emission； - GPU 加速； - 通过 OO/class hierarchy 构建复杂抽象层； - 将恒星背景预烘焙为 RGB 天球纹理； - 依赖单个固定坐标位置的“相机”。 可以保留 frozen-snapshot 模式作为 nmesh 内部的廉价在线 diagnostic，但它与本 repo 的最终 4D renderer 是不同层次的工具。 --- ## 3. 语言与总体实现风格 首选 **C**。 C++ 并不是当前项目的必要条件。需要的抽象主要可以通过： - `struct` - 函数指针表 - 显式 ownership - 大块连续数组 - thread-local workspace 完成。 性能敏感部分尽量采用： - bulk array processing； - SoA（structure of arrays）； - 避免 per-ray malloc； - 避免大量虚调用/函数指针 dispatch； - 避免细粒度锁； - 避免每个 triangle/ray 建 task； - 尽量顺序读取 metric 数据； - OpenMP 做 CPU 并行。 核心实现原则： > **bulk arrays + time-slab streaming + coarse-grained parallelism** --- # 4. 总体物理与数据流 整个程序不是“逐帧、逐像素独立 trace”，而是： 1. 预先生成相机 worldline 与 tetrad 轨迹； 2. 建立所有视频帧的初始 image-plane adaptive mesh； 3. 从所有帧收集当前 refinement level 需要的新 ray samples； 4. 将这些 rays 组成一个全局 `RayPool`； 5. 从视频结束时刻向过去，按 time slab 顺序加载数值时空； 6. 在每个 slab 内，把所有 active rays 一起推进到 slab 左边界； 7. ray 若在某个渐近端逃逸、达到红移暗阈值、耗尽预算或遇到数据/积分失败，则按类别终止； 8. 一轮 ray tracing 完成后，把 endpoint 数据回填到各帧 image mesh； 9. 根据局部 lens mapping 误差判断哪些 image-plane triangles 需要进一步细分； 10. 生成下一批新增 rays； 11. 重复若干 refinement passes； 12. 最终用局部 inverse lens maps 查询恒星 catalog； 13. 计算每个恒星像的位置、放大率和频移； 14. 将对应黑体 PSF 加到 HDR framebuffer； 15. tone mapping / 编码输出视频。 总结构： ```text Observer Track z^μ(t), e_(a)^μ(t) │ ▼ Movie / Frame Set │ ▼ Adaptive image-plane meshes │ ▼ Refinement analyzer │ ▼ new ray samples │ ▼ RayPool │ ▼ Spacetime time-slab stream (newest → oldest) │ ▼ ray endpoint: n∞, frequency shift, end/outcome/reason │ ├──────────► next refinement pass │ ▼ Final local lens maps │ ▼ Star catalog (n∞, T, amplitude) │ ▼ inverse mapping + μ + g + PSF │ ▼ HDR video frames ``` --- # 5. 为什么必须按 time slab 组织，而不是按帧组织 最终 BBH 4D 数据很可能超过单机内存。 因此不能对每帧独立： ```text for frame: load all required 4D metric data trace all rays ``` 否则同一段时空会被反复从硬盘读取。 更合理的是： ```text for time_slab = latest -> earliest: load metric data for this slab activate all rays whose camera time enters this slab evolve every active ray backward through this slab remove terminated rays discard slab ``` 这样： - metric I/O 主要取决于 NR 数据量； - 不随 4K 分辨率、帧数、supersampling 线性重复； - 数据访问可以是大块顺序读； - 适合本地 NVMe； - 适合未来超过 128 GB 的数据集。 --- # 6. Ray 生命周期与跨帧整合 不应一开始生成整个视频所有 rays 并常驻内存。 随着倒序时间扫描，在经过某个 frame 的 `t_camera` 时才激活该帧当前 refinement pass 的新 rays。 因此 active-ray pool 大小由以下因素决定： ```text 单位相机时间产生的 ray 数 × ray 在 NR domain 内的典型 coordinate-time 停留长度 ``` 而不是整个视频总 ray 数。 对于长时间停留在 strong-field region 的 rays，可以一直跨多个 slabs 保留。 --- # 7. Adaptive image-plane mesh 最终点源渲染不适合简单的“每 pixel backward lookup sky color”。 恒星在无穷远天球上是点源： \[ I(\hat n) \sim \sum_s F_s \delta(\hat n-\hat n_s) \] 因此需要先通过 backward ray tracing 得到局部 lens mapping： \[ F:\ (x,y)_\text{image} \to \hat n_\infty \] 然后在每个局部可逆 patch 上建立 inverse： \[ F^{-1}:\ \hat n_\infty \to (x,y)_\text{image} \] ### 基本单元：triangle 每个 image-plane triangle 的三个顶点都保存： - image-plane 坐标 `(x,y)`； - ray 的终态类别 `ESCAPED`/`DARK`/`UNRESOLVED`/`INCOMPLETE` 及其 reason； - 若 escaped：无穷远方向 `n_inf` 与所属 end； - frequency shift / redshift accumulator。 示意： ```c typedef struct { double x, y; double n_inf[3]; double log_g; uint8_t outcome; /* ESCAPED / DARK / UNRESOLVED / INCOMPLETE */ uint8_t reason; /* REDSHIFT_LIMIT / BUDGET_EXHAUSTED / ... */ uint32_t end_id; } LensVertex; typedef struct { uint32_t v[3]; uint8_t level; uint8_t flags; } LensTriangle; ``` --- # 8. 局部 adaptive refinement 局部 inverse map 的球面面积权重必须满足凸组合约束。当前实现以无符号 子面积除以总面积得到权重，在插值位置和频移前检查权重和有限、为正，且 `abs(sum(weights) - 1) <= 1e-8`；不满足时只跳过当前恒星在当前 triangle 中的像。通过检查后，将权重除以其总和，消除小的归一化误差。整 tile 包含的 catalog 快速路径也必须执行这项检查。 这是反演有效性保护，不是 refinement 阈值：细长 source triangle 的带容差 边测试可能误接纳外部点，无符号子面积之和便大于总面积，直接插值会把像 放到 image triangle 外。不能先归一化明显无效的权重来掩盖误接纳。 固定的无量纲容差 `1e-8` 依据 Schwarzschild 3840×2160、45° FOV、R=100、 coarse cell=32、level=4、J=0.2 网格上的诊断选定：正对中心及偏离 0.00161°、方位为 1°/44°/23° 的单星测试中，严格包含样本的 double 权重和误差不超过约 `2.5e-13`，正对中心的越界样本最小误差约 `5.9e-4`。 这不保证临界曲线附近的离散映射已经收敛，也不保证数值亮环无缺口。 每个 triangle 只关心自身局部映射。 不需要先把整个像平面映射到天球、再全局分类“第几阶像”。 对一个 triangle： 1. 已有三个 corner rays； 2. 计算或已有 edge midpoint / center rays； 3. 比较真实 mapping 与低阶插值； 4. 检查 orientation consistency； 5. 检查 Jacobian 是否接近奇异； 6. 若误差过大则 subdivide。 这意味着： - Schwarzschild 下无需显式按 winding number 分类； - BBH 下也不需要定义全局“第几阶像”； - 每个局部 patch 只需保持局部单值可逆； - 多像会自然表现为同一个 source direction 被多个 image patches 覆盖。 --- # 9. 为什么 refinement 需要多 pass adaptive refinement 与 out-of-core time streaming 存在天然冲突。 只有当一条 ray 已经完整追到无穷远后，才知道某个 triangle 是否需要细分。 但此时早期 metric slab 已经释放。 因此不能在单次 time sweep 中随时创建新 ray 并从相机时刻重新追。 采用： ```text Pass 0: trace coarse samples through whole spacetime analyze maps Pass 1: trace only newly requested midpoint/center rays analyze maps Pass 2: trace only newly requested rays ... until converged ``` 每个 pass： - metric data 仍然只顺序扫描； - 只追踪新增样本； - refinement 通常若干轮即可完成。 --- # 10. 相机不是固定坐标点 相机不应由： ```c double x[3]; ``` 定义。 否则结果强依赖某个 NR 坐标系中的“固定点”。 应由独立 trajectory generator 预先生成： \[ z^\mu(\tau), \qquad e_{(a)}{}^\mu(\tau) \] 其中： - \(z^\mu\)：observer worldline； - \(e_{(0)}{}^\mu=u^\mu\)：observer 四速度； - \(e_{(i)}{}^\mu\)：相机空间 tetrad。 trajectory generator 可以采用不同物理规则： - timelike geodesic； - Fermi-Walker transport； - parallel transport； - 人为给定加速度； - 人为叠加 camera pointing/rotation。 电影 renderer 本身只读取轨迹；单张入口可以用下述参数直接构造一个事件处的完整 tetrad。 示意： ```c typedef struct { double t; double tau; double x[3]; double tetrad[4][4]; } ObserverSample; typedef struct { ObserverSample *samples; size_t n; } ObserverTrack; ``` 每帧根据 `t_camera` 插值得到完整 observer state。 ## 10.1 单张的瞬时相机 单张与电影共享 `ObserverState` 和 ray 初始化。单张不是只指定三维位置： 由事件、坐标速度、指向与 roll 生成完整四速度和 tetrad。 `--observer-time T` 指定单张相机事件的坐标时（默认 $t=0$，接受任意有限值， 包括负值）；metric 求值、ray 初始化与 lens-map 帧元数据使用同一事件时刻。 此构造不要求静态时空或静止观测者；瞬时相机的 proper time 仍以零为参考。 `--observer-position X Y Z` 与 `--look-ra-deg` / `--look-dec-deg` 独立。 只有位置时取朝原点的坐标方向；只有指向时令 $\mathbf x=-R\mathbf d$。 `--observer-radius R` 仅用于缺省位置的推导（默认 30），与显式位置互斥。 几何参数全部省略时保留 Minkowski 原点、Schwarzschild 半径 30 的默认相机。 原点无法推导指向，要求显式 look；由位置推导方向且在极轴上时固定 RA=0。 RA 或 Dec 只给一个时，另一个取默认值（90°、−90°）。 `--observer-velocity VX VY VZ` 表示 $v^i=dx^i/dt$，默认零。 由事件处的 $\alpha,\beta^i,\gamma_{ij}$ 计算 $$ Q=-\alpha^2+\gamma_{ij}(v^i+\beta^i)(v^j+\beta^j),\qquad u^\mu=\frac{(1,v^i)}{\sqrt{-Q}}. $$ 要求 $Q<0$，未来指向选择 $u^t>0$。不能用欧氏 $|\mathbf v|<1$ 或 $r>2M$ 代替类时性检查；删除旧的 `--observer-inward-speed` 参数。 构造失败直接报错，不修改速度或退回静止相机。 采用标准 ICRS 坐标种子 $$ \mathbf d=(\cos\delta\cos\alpha_{\rm RA},\cos\delta\sin\alpha_{\rm RA},\sin\delta), $$ 将 $D^\mu=(0,\mathbf d)$ 投影为 $F^\mu=D^\mu+u^\mu(u_\nu D^\nu)$，归一化得到 forward。 随后把天球北向和西向的坐标种子投影、正交化为 up、right，保留标准画面定向。 因此 look 表示坐标方向在相机静止空间内的投影，不保证特定远方物体位于画面中心。 `--camera-roll-deg` 默认零，正角按 $e'_2=\cos\rho\,e_2+\sin\rho\,e_3$、 $e'_3=-\sin\rho\,e_2+\cos\rho\,e_3$ 定义。 observer 构造只接收当地 metric 和已补全的参数，不加载 slab、不分类 ray。 调用方在昂贵的 catalog/PSF 初始化前验证相机及 backend 数据域。 相机合法性只由 metric 可用性、四速度 timelike、时间定向和 tetrad 正交归一决定； 视界内、旧 cutoff 内的相机都是正常渲染目标，位置本身不决定 ray 终态。 向过去追踪的高红移阈值截断仍正常生效（见 §18）。 单张相机参数与轨迹输入、lens-map 导入互斥；导入仍跳过 metric 与 observer 初始化。 验证包括 tetrad 正交归一和 null 初始化、平直时空平移不变性与解析光行差/多普勒、 KS 视界处和视界内的向过去径向逃逸及频移，以及相同 observer 的单张/电影单帧一致性。 --- # 11. Ray 初始化 给定相机 film coordinate `(x,y)`： 1. camera projection 得到 observer local frame 中的 photon direction； 2. 与 observer tetrad 相乘得到 spacetime 中的 \(k^\mu\)； 3. 转换到 geodesic integrator 所使用的 3+1 ray variables。 例如： \[ k^\mu = e_{(a)}{}^\mu k^{(a)} \] film model 只负责： \[ (x,y) \to \hat n^{(i)}_\text{camera} \] 它与 gauge、metric backend、observer worldline 都解耦。 --- # 12. Ray 状态 采用 coordinate time \(t\) 作为 ODE 自变量，以便天然和 time slabs 对齐。 单 ray 概念结构： ```c typedef struct { double t; double x[3]; double Pi[3]; double log_E; double h; uint32_t sample_id; uint32_t frame_id; uint8_t status; } Ray; ``` 真正批量计算时建议 SoA： ```c typedef struct { double *restrict t; double *restrict x0; double *restrict x1; double *restrict x2; double *restrict p0; double *restrict p1; double *restrict p2; double *restrict log_E; double *restrict h; uint32_t *sample_id; uint32_t *frame_id; uint8_t *status; size_t n; size_t capacity; } RayPool; ``` 目的： - cache locality； - SIMD； - OpenMP； - stream compaction； - 避免 per-ray allocation。 --- # 13. Spacetime backend 所有上层渲染逻辑都不应该知道 metric 来自哪里。 需要支持至少： - Minkowski； - analytic Schwarzschild； - 未来可加 analytic Kerr； - numerical nmesh 4D spacetime。 核心不是“单点 metric 对象”，而是**当前已加载 time slab 的 evaluator**。 概念接口： ```c typedef struct SpacetimeSource SpacetimeSource; typedef struct MetricSlab MetricSlab; typedef struct { int (*load_slab)( SpacetimeSource *src, double t_hi, size_t memory_budget, MetricSlab **out); void (*free_slab)(MetricSlab *slab); int (*eval)( const MetricSlab *slab, double t, const double x[3], void *workspace, MetricData *out); int (*classify)( const MetricSlab *slab, double t, const double x[3]); } SpacetimeOps; ``` 解析 backend 的 `load_slab()` 可以近似 no-op。 nmesh backend 则真正加载一段时间的 3D outputs。 --- # 14. 统一的 3+1 metric 数据 初步采用： \[ \alpha,\quad \beta^i,\quad \gamma_{ij},\quad K_{ij} \] 作为主要 3+1 变量。 空间导数： \[ \partial_i \alpha,\quad \partial_i\beta^j,\quad \partial_i\gamma_{jk} \] 由 backend evaluator 提供。 示意： ```c typedef struct { double alpha; double beta[3]; double gamma[6]; double K[6]; double d_alpha[3]; double d_beta[3][3]; double d_gamma[3][6]; } MetricData; ``` --- # 15. 时间导数 时间导数不能简单假设全部从文件直接输出。 应按量的性质选择来源。 例如： \[ \partial_t \gamma_{ij} = -2\alpha K_{ij} + \mathcal L_\beta \gamma_{ij} \] 可以直接由 3+1 内部量构造。 对于 lapse/shift： - 若 gauge evolution equations 已知，可从 gauge RHS 构造； - 或从 temporal interpolant 解析求导。 原则： > 若已经建立 \(u(t,\mathbf x)\) 的时间插值，就从同一个 interpolation polynomial 解析求 \(\partial_tu\)，而不是另做 finite difference。 且不要每个 geodesic RHS 都计算用不到的时间导数。 --- # 16. nmesh 数值时空 backend nmesh 当前 3D 输出格式已经非常适合 renderer： 对每个需要输出的 timestep： - 每个 AMR leaf node 一个文件； - 文件名按 node 的树状 refinement 命名； - 从文件名可恢复 refinement tree； - 文件内保存每个 nodal point 的 global coordinates； - 保存所选物理量。 nmesh： - DG； - node 内用 local coordinates 上 Legendre nodes 对应的 Lagrange basis； - production 中每边一般最多约 12 个 nodal points； - AMR 为 h-refinement； - 典型网格为中心六面体 + 外层 cubic-sphere shells。 因此读入一个 node 文件后，可以直接构造 node 内： \[ u(\xi,\eta,\zeta) = \sum_{ijk} u_{ijk} L_i(\xi)L_j(\eta)L_k(\zeta) \] 空间导数直接对 Lagrange polynomial 求导。 不需要： - resample 到 Cartesian uniform grid； - 额外输出 derivative fields。 --- # 17. nmesh time interpolation h-AMR 可能使相邻 timestep 的树结构不同。 但不要求时间方向上 node 一一对应。 对任意固定 global point \(\mathbf x\)： \[ u_n(\mathbf x) = I_{\mathcal M_n}[u](\mathbf x) \] \[ u_{n+1}(\mathbf x) = I_{\mathcal M_{n+1}}[u](\mathbf x) \] 先分别在各自空间 mesh 上求值，再做： \[ u(t,\mathbf x) = T[u_n(\mathbf x),u_{n+1}(\mathbf x),\ldots] \] 因此： - 每个 time slice 可独立恢复空间函数； - AMR topology 不需要跨时间匹配； - time slab 只负责保存足够多邻近 slices 供时间插值。 slab 边界需要少量 overlapping temporal ghost slices。 --- # 18. 黑洞终止与暗终态 目标使用 moving-puncture BBH，而不是 excision。正常终态**不使用** horizon 内位置 cutoff、AH-calibrated puncture 小球或 armed/re-entry 状态机判定物理捕获。过去向光线 围绕渐近端逃逸、能量阈值截断和经可靠识别的渐近轨道组织；达到红移阈值后停止、渲染 为黑，是已确定需求。 形式（相机相对局域能量增长，对全部 backend 统一；具体阈值通过小型 oracle/convergence test 标定，不宣称由论文给定）： \[ L-L_0=\ln\!\frac{\alpha p^0}{(\alpha p^0)_0}\ge L_\text{dark} \Rightarrow \text{DARK (redshift limit)} \] `L_0` 是相机事件的参考值，对全部 spacetime backend 生效，并随 ray 状态跨 slab 与 retry 传递；减去 `L_0` 只改变判据参考，不重置光子能量或频率比 `g`。必须区分 `L` （Eulerian 能量对数）、`ln(p^0)=L-ln(alpha)` 和真正连接源端得到的 `g=E_camera/E_source`。当前默认 `L_dark=8`，由 CLI 参数覆盖，不宣称为论文值。 终态分为四类（详见 §18A 与终点协议）： - `ESCAPED`：成功完成某个 end 的外推，必须带有效 `end_id`、`n_infinity`、`g`； - `DARK`：无天空贡献的暗终态，当前主要为红移阈值截断 `REDSHIFT_LIMIT`；不同 dark reason 不制造 mesh seam； - `UNRESOLVED/BUDGET_EXHAUSTED`：轨迹仍可信但计算预算用尽，可重试； - `INCOMPLETE/FAILED`：历史耗尽、域外、invalid metric、I/O、积分误差不可控、 unsupported chart 或 protocol error；不得伪装成 dark。 数值失败、单次 metric eval failure 或单个超阈值 trial step 均不得推断为物理 capture； 阈值只检查可信的初始或 accepted 状态。`ASYMPTOTIC_TRAPPED` 与 `SINGULARITY` 仅在存在 可靠 backend 判据及明确源边界条件时启用。有限几何分辨率导致的 shadow 略偏大与未解析 高阶像尾部由三角形近似处理，并保留 triangle provenance 与面积统计；这属于渲染近似， 不是物理捕获。跨 chart、跨 region 或穿越视界本身不是暗终态；moving-puncture trumpet 不解释成可穿越的第二个宇宙。production renderer 不假定每次 NR run 都运行昂贵的 AH finder，也不依赖 capture sidecar。 解析 Schwarzschild 的 policy version 3 使用相机相对局域能量增长 \[ A_0 = L - L_0 \ge L_\text{dark}. \] `L_0` 是该 ray 积分起点的参考值，随 ray 状态跨 slab 与 retry 携带，不在每次 检查时用即时状态重算，也不重置光子能量或 frequency ratio。Eulerian 观者测得 的能量为 `e^L`，故 `A_0 = ln(E_euler/E_euler,0)`；一个常数相机 boost 在减法中 抵消，因此大 boost 或内部相机不会仅因初始 `L` 大而被判暗。这正是与视界无限 红移对应的局域相对量。对静止时空，理论上也可用 Killing 相对量 `A_K = L - ln|E_K|`，但沿精确光线 `A_K - A_0 = -ln|alpha_0 - beta_0.Pi_0|` 只是初始常数，不能仅凭守恒证明其优于 `A_0`；该 stationary 参考作为独立对照保留 在 `test_termination_oracle.c`，不进入生产判据，也不移植到动态 NR。阈值 `L_dark=8` 仍是待标定参数；生产条件不含绝对 `L` 或 backend applicability guard。 所有 midpoint probes 同样属于完成性检查范围。失败/未决 probe 保存为 off-mesh witness；未决 witness 合并到下一轮续追请求，完成 witness 可复用。重试回填使相关 叶子决策失效，不能借用其他叶子的 `evaluated` 标记跳过新出现的边界。最大层数、 最长边和面积共用同一停止判据；近似标黑统计报告实际最大边长、最大面积及层数停止 数量。Replay 使用文件保存的几何策略进行同样的完成性检查，诊断覆盖必须报告 `INCOMPLETE`，不得绕过发布 gate。conformity/几何限制取消全部请求边的 triangle 必须显式 settle（记录 `evaluated`），不得每代重复请求同一组 probes；被取消而 几何仍允许的 UUD/UDD 结算为 budget-incomplete，达到停止尺度的结算为近似标黑。 witness 提升为正式 midpoint 时原地复用同一 vertex id，只保留一份连续状态；只有 真正 off-mesh 的未决 witness 独立重试。lens-map 每帧保存累计 `retry_requests`， provenance 保存 coordinate-time step 与初始 step 预算，使实际积分来源与成本可 replay。 诊断输出所有权：所有 ray 失败报告只在 CLI/main 的串行后处理阶段发出， ray-tracing 的 OpenMP worker 与 geodesic 热循环不写任何日志（既有 PSF splat 诊断不属于此范围）。稳定的 `ray_reason_name` 覆盖全部 `RayReason`；默认（非 verbose）也按 frame/phase 输出 `INCOMPLETE` 的 reason 直方图与真实 frame id、相机坐标时间，verbose 或 GR_DEBUG 才追加每 reason 有界 数量的代表性样本（generation/phase、sample id/kind、持久 vertex id、film 坐标、 可信 endpoint 的 stop 状态与 accepted/rejected/RHS 计数）并报告被抑制数量。 `UNRESOLVED/BUDGET_EXHAUSTED` 只在帧边界统计确认存在 blocking triangle 时按帧 报告，与数值 `INCOMPLETE` 区分；replay 使用文件中保存的几何与顶点诊断，不虚构 stop 状态，也不把未重试的旧失败当作新失败重复计数。 `RayReason` 的诊断粒度是追加式扩展：原 coarse reason 的数值 0..9 冻结不变， 更细的失败原因一律追加在 `RAY_REASON_IO_ERROR` 之后，`RAY_REASON_COUNT` 只作 计数/未知 sentinel 且不写入 wire。lens-map 的二进制布局（v2/v3）与字段偏移不 改变，新 reader 读旧值原样、读追加值按新名解释，并拒绝 `>= RAY_REASON_COUNT`； 旧 reader 因 coarse 校验仍会安全拒绝它不认识的新值。每个细原因通过纯函数 `ray_reason_category()` 归入既有 coarse 类别（PROTOCOL_ERROR、 INTEGRATION_ERROR、UNSUPPORTED、IO_ERROR 等），供只需要粗分类的调用方使用。 该扩展只细化失败状态的命名与上报，不改变任何物理终态判据、终止分类或 retry 规则。 --- # 18A. 渐近外区、escape worldtube 与 endpoint 协议 本节冻结“到达 escape 区域就终止”这一旧行为被替换后的职责边界，是该协议的权威 约定与唯一长期记录。 ## 18A.1 问题 旧判定只看光线当前位置是否在某个 escape 半径之外，不看传播方向。因此当相机 本身位于 escape 球外时，所有光线在初始化后立即被判为逃逸，连本应进入强场区的 光线也不会积分。旧实现还把“外推到无穷远”简化为“在 escape 球处删除”，频移因此 带有 $O(M/R)$ 的误差。 ## 18A.2 ray 生命周期三段 1. **相机位于 escape worldtube 外**：由公共渐近外区模块判断光线是否会与 worldtube 相交。 - 相交：把光线外推到第一次由外向内穿越，并从该 entry event 开始交给 backend 内区积分； - 不相交：直接把光线外推到对应无穷远天球，写 `ESCAPED` endpoint。 2. **光线在 backend 内区积分**：不再因为“当前位置处于 escape 区域”立即终止； 只有沿 ray 的过去传播方向发生**有方向的 inside -> outside 穿越**时才进入外区 收尾。 3. **穿越后**：公共渐近外区模块把有限半径处的 canonical photon state 推到无穷远， 得到 `n_infinity` 和 `frequency_ratio`。 一个 backend 可以声明多个渐近远端；当前实现只暴露一个 `end_id`，但 endpoint 与 ray 状态中不得把“整个时空只有一个无穷远”写死。 ## 18A.3 职责边界 backend 负责声明： - `end_id` 及其稳定编号； - 外区模型种类：`MINKOWSKI` 或 `SCHWARZSCHILD_MONOPOLE`； - 从该远端看见的渐近质量 `mass`（允许为零）； - 渐近参考系的 origin 与空间基（在 backend 坐标中表达）； - 给定 coordinate time 的 escape worldtube 球心 `center`、速度 `velocity`、 半径 `radius`、半径变化率 `radius_rate`； - worldtube 描述有效的时间区间与运动分段边界； - backend 内一点属于哪个候选 end 的 outer region，或当前无法分类。 backend **不**负责：球外传播、entry/miss 判定、无穷远方向、pre-route 或 endpoint 写回。 公共渐近模块负责： - 将 backend photon state 与统一 canonical state 双向转换； - 球外相机的 entry/miss 判定； - 对 miss 光线直接生成 infinity endpoint； - 把 entry 光线传播到 worldtube 的第一次由外向内穿越； - 把内区积分产生的由内向外穿越传播到无穷远； - $M=0$ 使用精确 Minkowski 几何；$M>0$ 固定同心球使用内建 monopole 近似； - 返回明确状态码，而不是用 NaN 或任意 fallback 掩盖适用域错误。 geodesic/ray 生命周期层负责： - 初始化时调用 pre-route； - 保存 entry event，并在 slab sweep 到达 entry time 时激活内区积分； - 每个 accepted ODE step 后检测有方向的 crossing 并局部化第一次根； - 调用公共外区模块完成 endpoint。 ## 18A.4 canonical photon state 与时间约定 canonical state 至少包含：coordinate time $t$、渐近参考系中的位置、传播方向与 能量/频移所需的 photon momentum 信息、以及 `end_id`。 renderer 沿过去方向积分。文档中使用的空间单位方向唯一约定为**过去传播方向** $\mathbf w$：令 $s=t_\text{camera}-t\ge 0$，则局部轨迹满足 $\mathbf x(s)=\mathbf x_0+s\,\mathbf w+\dots$。`n_infinity` 是光线在无穷远处的 来向，即 $\mathbf w$ 在 $s\to\infty$ 的极限，因此与旧实现中 `normalize(-gamma^{ij} Pi_j)` 的符号约定一致。 在静止时空（Minkowski 与 Kerr–Schild Schwarzschild）中，沿测地线守恒的 photon 能量为 \[ E_\infty=-p_t=\alpha p^0\left(\alpha-\beta^i\Pi_i\right), \] 其中 $p_i$（即 `Pi` 的协变版本）满足 $p_i=\alpha p^0\,\Pi_i$。相机归一化取 $E_\text{camera}=1$，故 \[ g=\frac{E_\text{camera}}{E_\infty}=\frac{1}{\alpha p^0(\alpha-\beta^i\Pi_i)}. \] 旧实现对 $M>0$ 在 escape 球处直接返回 $\exp(-\log(\alpha p^0))$，是上式在 $\beta^i\Pi_i\to0,\alpha\to1$ 下的近似。 ## 18A.5 worldtube 与穿越方向 球面 worldtube： \[ F(t,\mathbf x)=|\mathbf x-\mathbf c(t)|^2-R(t)^2. \] $F>0$ 外、$F<0$ 内、$F=0$ 边界。边界点（$F=0$）必须结合过去传播方向的斜率 $\mathrm dF/\mathrm ds$ 分类：$\mathrm dF/\mathrm ds<0$ 视为即将进入、$\ge0$ 视为 向外或切触；因此相机恰在 $F=0$ 且 past-inward 才按 INSIDE 处理，past-outward 与 tangent 都按外层 route 处理。inside -> outside crossing 要求 $\mathrm dF/ \mathrm ds>0$ 的严格符号变化（$F_{\rm before}\le0$ 且 $F_{\rm after}>0$）； 仅有 $F_{\rm after}=0$ 的单点切触不算 crossing，需等下一步是否真正到 $F>0$。 必须区分方向： - camera pre-route 的 entry 是沿过去传播方向第一次 outside -> inside； - 内区 escape 是沿过去传播方向第一次 inside -> outside； - 某次采样发现 $F\ge0$ 不能独立构成 escape。 对步进端点接近零、切触和跨越 motion-segment 边界，使用显式容差和有界 root localization；不得用固定位置 epsilon 把 tangent 误判成 crossing。 ## 18A.6 $M=0$ 外区 渐近惯性系中为解析直线传播。固定球用 ray-sphere 二次方程取沿过去传播方向最早 的合法根；匀速移动球在分段内把球心写成 $\mathbf c(t)=\mathbf c(t_0)+\mathbf v(t-t_0)$， 令 $\mathbf d=\mathbf x_0-\mathbf c(t_0)$、$\mathbf q=\mathbf w+\mathbf v$，entry 满足 $|\mathbf d+s\mathbf q|^2=R^2$（半径线性变化时右端为 $(R_0-R_\text{rate}s)^2$）。 任意加速球的未来接口使用分段 bracketed root driver；若 backend 历史在判定完成前 结束，返回 `TIME_RANGE_EXHAUSTED`，不得武断判为 miss。 $M=0$ 的 finish 是平凡的：$\mathbf n_\infty=\mathbf w$、 $g=\exp(-\log(\alpha p^0))$（flat 中守恒）。 ## 18A.7 $M>0$ Schwarzschild-like 外区 首版严格限制：球心固定、escape 球与 monopole 同心、$R/M\ge64$、$M>0$、相机与 worldtube 位于该外区。不满足则返回明确的 unsupported/domain 状态；不静默退回 Minkowski，也不把一般移动 Schwarzschild 球解释成瞬时静态球。 无量纲量 $\rho=r/M$、$\beta=b/M$， \[ Q(\rho,\beta)=1-\beta^2\frac{1-2/\rho}{\rho^2}. \] escape 球处切触阈值 $\beta_R=\rho_R/\sqrt{1-2/\rho_R}$。entry/miss 的拓扑分类优先 使用解析阈值和方向信息，不由低精度查表决定。 角度 primitive 为过去传播方向从半径 $\rho$ 到无穷远扫过的单调外向方位角 \[ \Phi(\rho,\beta)=\int_\rho^\infty \frac{\beta}{\rho'^2\sqrt{Q(\rho',\beta)}}\,\mathrm d\rho' =\int_0^{1/\rho}\frac{\beta\,\mathrm du}{\sqrt{1-\beta^2u^2+2\beta^2u^3}}. \] turning radius 满足 $\beta^2=\rho_\text{turn}^3/(\rho_\text{turn}-2)$。守恒的 impact parameter 与角动量满足 \[ \beta=\frac{|x\times\Pi|}{\alpha-\beta^i\Pi_i},\qquad \mathbf N=\widehat{x\times\Pi}, \] 无穷远方向由 $\hat{\mathbf r}=x/|x|$ 绕 $\mathbf N$ 旋转 $\Phi(\rho,\beta)$ 得到。 turning map 与 Schwarzschild coordinate-time transfer 使用离线验证过的有界 residual、Chebyshev 表或解析主项；运行期不得建表。 若 time-transfer 无法在声明域内满足误差标准，则保留 $M=0$ 实现和接口，不把未经 验证的时间公式写入生产代码，也不得降低验收标准。 ## 18A.8 nmesh outer-shell 约定（仅约定，不实现） - 最外层必须是有明确六个面的 cubed-sphere shell； - outer boundary 在渐近 frame 中是固定中心、固定半径球面； - 提供 $R$、对应 end 的质量和 frame metadata； - Schwarzschild monopole 模式要求 $R/M\ge64$；$64M$ 只是拒绝更靠内 junction 的硬下限。nmesh 有 AMR，生产数据应把 outer shell 放到尽可能大的半径，建议以 至少接近解析 backend 当前的 $256M$ 为目标； - DG element 边界本来允许场跳变，故 worldtube junction 不要求两侧 metric pointwise 连续； - 穿越时匹配 boundary local tetrad 中的 photon direction/energy，并用分辨率与 outer-radius convergence test 验证，而不是强行匹配坐标分量。 ## 18A.9 失败语义 遇到以下情况必须返回明确状态并上报，不得 fallback： - worldtube 历史不足，无法判断 first entry； - Schwarzschild 外区不是固定同心球； - $R/M<64$； - canonical state 无法保持 null constraint 或 round-trip 精度； - time coordinate 约定不明确； - 高精度 reference evaluator 在域内不收敛； - 表在 seam 或 grazing 区域超过误差限。 ## 18A.10 Backend 能力与路由约束 - 公共接口以 end descriptor、escape worldtube sample、canonical photon state 与 endpoint `end_id` 分离远端声明、几何事件、状态转接与终态来源； - $M=0$ 固定球与匀速移动球的 camera pre-route、directed inside -> outside crossing 与 `PENDING_ENTRY` 生命周期采用同一几何约定，适用于 Minkowski 与 Alcubierre 的平直外区； - 单帧与 movie 必须共用同一 pre-route/生命周期实现； - $M>0$ 固定同心 Schwarzschild monopole 外区以**解析 Carlson 椭圆积分**实现 （见 18A.11），接入解析 Kerr–Schild backend；相机位于 escape 球外的 pre-route、entry 传播、directed crossing 与 infinity endpoint 均可用； - worldtube 的运动能力必须显式声明。分段匀速、线性半径变化可在各段内求首次 crossing，并在段边界重新采样；任意加速运动在没有可信区间界时必须报告 `UNSUPPORTED`，不得把离散采样没有发现根解释为正常 miss。扩展加速运动能力需要 同时提供可信的段内界或专用首次 crossing 解法。 - `end_id` 贯通 endpoint、`LensVertex`、`terminal_mismatch` 与 `discrete_jacobian`，避免跨 end 插值；lens-map 保存 `end_id`。不同物理 region 的 reachable-end 路由及 catalog 绑定仍需独立定义， 不能把能够保存 end 标识解释为已支持完整解析延拓。 - 生命周期为三态：`NO_ENDS`（才允许 legacy）、`DIRECTED_READY`、 `PROTOCOL_ERROR`（descriptor 读取失败或不支持的 exterior，显式 `INVALID/UNSUPPORTED`，绝不退回 legacy）。pre-route 也先校验所有 descriptor 与 exterior kind，再判定 worldtube 内外，因此“声明了不受支持 exterior 而相机 恰在球内”不会被静默接受。 - 内区积分时若 worldtube sample 变为 `valid == 0`（历史耗尽），作为一等 terminal reason 返回 `INCOMPLETE/TIME_RANGE_EXHAUSTED`、保留 `end_id`，并由 `RayPool` 记为 `RAY_POOL_TERMINATED`，与 pre-route 的耗尽语义一致。内区 crossing localizer 定位过程中任何 sample 失败都直接传播具体状态，不返回一个正常 entry。 - Minkowski entry quadratic 用稳定根公式（$q=-\tfrac12(b+\mathrm{copysign} (\sqrt\Delta,b))$，取最小正根）；$c=0$（相机在边界）时按 $b=\mathrm dF/ \mathrm ds$ 分类。worldtube sample 必须有限且 $R>0$，否则 `INVALID`。 Alcubierre 等平直外区共用此解法；相对位置/速度、系数、判别式和根采用 `long double` 中间量，最后才转换为 backend 的 double 状态。某些平台的 `long double` 不提供额外有效位，因此精度提升不能替代入口验证。 送入根公式前，以系数最大绝对值的二进制指数共同缩放 $a,b,c$（最大值进入 $[1/2,1)$），避免判别式乘积溢出；不改变传播参数及根，也不覆盖系数构造前 已发生的误差。非有限系数或缩放后非零系数落入 subnormal 范围，转入不确定 路径而非当作退化线性式。判别式 $b^2-4ac$ 和入口斜率采用普通 `long double` 运算，不使用显式 FMA；精度/成本对照见 `benchmarks/quadratic_precision/`。 本机软件 `fmal` 成本明显，测试未显示足以抵偿成本的精度收益；不确定带、 几何验证及后备定位仍保留，不能把此样本结论当作全参数可靠性保证。 - **入口快路径与后备定位分离**：解析相交只提供候选入口；采用前须在实际 backend 坐标、实际入口时间重新检查 worldtube 残差。通过原有几何舍入容差的候选沿快路径 激活；未通过的候选沿同一外区 geodesic，从原相机事件重新求值并用确定在外/ 严格在内的区间定位首次入口。公共 localizer 不依赖具体 metric backend， 各已支持的外区提供轨迹求值与首次入口 bracket；不能把离散采样无交点当作 miss。 miss 判定必须排除真实入口；擦边或求根病态造成的数值不确定性不能作为 `ESCAPED` 的依据，允许保守地生成候选入口并进行后备验证。 浮点判别式等于零不构成精确擦边证明；重建最小点的单次正残差也不构成 miss 证明（大时间/坐标的舍入会移动该最小点）。精确擦边可由独立几何证据 排除入口，例如固定半径段内某个不变坐标的精确距离已不小于球半径； 无法获得这种证据且无法构造可信 bracket 时仍返回 `ENTRY_UNCONFIRMED`。 后备路径返回 bracket 收敛后的可表示内侧轨迹状态，不做径向位置投影，也不放宽 原 worldtube 检查容差。未能确认首次入口须报告 `INCOMPLETE/ENTRY_UNCONFIRMED`， callback、历史及几何错误仍传播各自具体 reason，不能改写为 capture 或 escape。 定位在 pre-route 中完成，不在 slab sweep 中倒回相机时间；不重置相机 `L0`， 内区 accepted-step/lookback 预算仍从最终激活事件起算。成功的后备入口由 `entry_fallback_evaluations` 记录 localizer 的轨迹求值次数（不含构造 bracket 的探测），不混入内区 RHS 计数；临时状态由 调用者独占，不引入共享可变缓存。此标记不改变 lens-map 二进制布局。 - **构造期验证优先**：`spacetime_create_*()` 成功即承诺该 source 已可安全光追。 每个 constructor 在安装 ops/context 后调用公共 `spacetime_source_finalize()`： 检查 ops/context 完整、`end_id` 唯一且非 `NONE`、exterior kind 受支持、 mass/frame 有限。**完整 worldtube 历史（所有 segment 边界、中心/半径连续性、 全程 $R>0$）由 backend constructor 负责**，普通解析 backend 在参数校验中完成； 失败时 constructor 销毁 context 并返回错误，不存在半构造可用的 source。 - **motion-segment 有效域**：候选二次根必须满足 $0\le s\le s_{\rm segment}$； 本段无有效根时推进到下一段边界（用 `nextafter(boundary,-∞)` 进入下一段）重新 求根；只有最后一个向过去开放的恒速段无根时才判 `ESCAPED`；segment 预算耗尽 属于内部失败，返回 `INVALID`。因为 constructor 已保证全程 $R>0$，正常 routing 不再做 segment-collapse 判定；但根处仍保留一次 $R_0-\dot R\,\sigma>0$ 检查（一次乘减），用于防御绕过 constructor 的 backend， 避免负半径伪 entry。 - **运行期防御**：所有 backend（含 Schwarzschild）都经公共 `worldtube_sample()` 读取，并保留一个便宜的 callback trust-boundary 检查 `isfinite(radius) && radius > 0`：callback 失败 → `INVALID`、`valid=0` → `TIME_RANGE_EXHAUSTED`、NaN/非正半径 → `INVALID`。这是防止第三方 backend 或 测试绕过 constructor 的安全网，不再承担正常配置验证；containment 循环原样 传播这些状态并保留 `end_id`。 - turning radius 用 bracketed bisection + Newton polishing。相机紧贴大球时的 entry 精度受 $|d|^2-R^2$ 输入条件数限制，由 entry-time 预算覆盖。 - generic accelerated driver 的初始条件同样用 $F<0$，或 $F=0$ 且 $\mathrm dF/\mathrm ds<0$ 才计为 entry；搜索只用严格 $F<0$ 的采样点确认 entry，单点 $F=0$ 切触不算。 - 终态契约已升级为两组枚举：渲染类别 `ESCAPED/DARK/UNRESOLVED/INCOMPLETE` 与诊断 reason。位置 cutoff 与 `SPACETIME_RAY_CAPTURED` 已删除；正常暗终态为 相机相对局域能量增长 `L - L0 >= L_dark`（默认 `L_dark=8`，可用 `--dark-threshold` 覆盖，对全部 spacetime backend 统一生效）。计算配额耗尽 返回可重试的 `UNRESOLVED/BUDGET_EXHAUSTED`， 保留最后可信连续状态；数据/积分/历史耗尽返回具体 `INCOMPLETE` reason。 - `eval`/`eval_slab` 返回 `SpacetimePointStatus`（时间不足、域外、invalid metric、内部错误）；observer 合法性只由 metric 可用性、timelike 四速度、 时间定向和 tetrad 正交归一决定，不再调用位置分类。 - 三角形决策实现 $E/D/U$ 账目：`UUU` 与含逃逸顶点的未决组合强制追加预算重试； `UUD/UDD` 在几何停止尺度处近似标黑并记录 image-plane 面积与 triangle provenance；含 `INCOMPLETE` 的三角形计为错误，不参与近似标黑。总资源上限 耗尽而未决时报告 incomplete，诊断可用 `--allow-incomplete` 覆盖。 - lens-map 文件格式升级为 v2：显式存储 `end_id/outcome/reason`、triangle `approx_black` 以及 policy/retry/几何阈值 provenance；v1 文件被明确拒绝。 ## 18A.11 $M>0$ 解析外区：闭式约化与验证 阶段 C 采用的不是 Chebyshev 表，而是把外区积分闭式约化到 Carlson 对称积分， 并在 double 下验证。记 $\rho=r/M$、$u=1/\rho$，$P(u)=1-\beta^2u^2+2\beta^2u^3$。 - **角度 primitive.** 对三次 $P$ 的实根分支用实 Legendre 形式 \[ \Phi=\frac{\sqrt2}{\sqrt{C-A}}\left[F(\varphi(u),\kappa) -F(\varphi(0),\kappa)\right],\qquad \sin^2\varphi=\frac{u-A}{B-A},\quad \kappa^2=\frac{B-A}{C-A}, \] $A max_step` 直接返回 protocol error，不偷偷填默认值。初始 trial step 复用 `coordinate_time_step`。 - DP54 采用标准 Dormand–Prince 5(4)7M tableau，接受 5 阶解，误差用 `h*sum((B5-B4)*k)` 直接计算；位置误差尺度 `atol_x + rtol*max(|dx|,|h*k1.x|)`，`Pi` 用 `atol_Pi + rtol*max(|Pi_before|,|Pi_candidate|)`， `L` 用 `atol_L + rtol*max(1,|dL|)`；控制器 safety `0.9`、factor clamp `[0.2,5]`、 exponent `1/5`，无 FSAL（每次 trial 7 次 RHS），拒步不提交状态。 - 自适应控制状态属于各条 ray：保存 `integration_start_time`（activation/entry 时间， 非相机能量参考 `L0`）、`next_step`、累计 rejected/RHS 计数与 `previous_rejected`； `steps` 表示 accepted 步数。状态须跨 slab 与预算重试保留；endpoint 保存最后可信 轨迹及控制/成本状态，使续追无需从相机重放，也不重置 `L0`。 - 控制/成本状态在 scheduler 各层按值传递：batch 路径的 `RayPool` SoA 保存同一组 per-ray 字段，slab activation 只重置新 ray 的 accepted 步数，绝不重置 adaptive control state；frame/movie 的续追请求以完整 resume state 重建（frame 与 movie 共用同一 helper，禁止逐字段手抄或漏字段）。endpoint 记录本次实际授予的 time budget，供 retry 层独立累计。off-mesh witness 原地提升必须整体复制该状态， 重试回填保留连续状态与成本计数。 - 计算配额与回溯时间是两条独立、各自饱和的 retry 预算：step 增量与 time 增量互不 替代，只有"实际耗尽且仍可增长"的 quota 才允许发起 retry；任一 blocking quota 已 达硬上限时不得仅凭另一 quota 产生无效重试，直接按 budget-incomplete 处理。两个 quota 的实际配置总 grant 必须与已 spent 数量分开记录，retry 采样取实际 grant， 不得由已花费步数反推授权。"可增长"指新 quota 必须打开严格更大的可表示区间：整数 步长严格增加，time 左边界 `start - quota` 严格前移；因大 time origin 或增量过小 而 round 到同一区间的增量不算有效增长，不得发起永远无法推进的 retry。两者都到上限 后 UUU 仍为 budget-incomplete，UUD/UDD 在几何停止尺度可按 finite-resolution 近似 标黑。自动补齐的 time 重试增量由 `max_lookback_time` 指定，总上限为其四倍 （对有限数值范围饱和）；这是可覆盖的 resource quota policy，不是物理终态条件。 - 每步对 slab 左边界与显式回溯预算左边界截断；DP54 不以 `coordinate_time_step * max_steps` 定义时间范围。回溯预算耗尽而轨迹仍可信 返回 `UNRESOLVED/BUDGET_EXHAUSTED`；真实的 slab data time failure 仍是 `INCOMPLETE/TIME_RANGE_EXHAUSTED`。`max_steps` 仍作为 accepted 步数额度。 各 trial 先确定可表示的目标时间，再用实际 `target-before.t` 推进状态、构造 dense interpolant 与提交时间；不得以舍入前的步长推进空间、舍入后的步长记时间。 worldtube 采样的坐标时间同样使用实际差值回推运动，不假设半步总是可表示。 - 初始 accepted 状态的 RHS/metric 失败直接报告具体 point reason，不缩步；后续 stage 的 `OUT_OF_DOMAIN/INVALID_METRIC` 可有界缩步重试；`TIME_UNAVAILABLE/INTERNAL_ERROR` 不靠无限缩步；最小步长或时间不可进导致失败时报告 `INCOMPLETE/INTEGRATION_ERROR`。 拒步上限耗尽时保留具体 stage 原因，或报告积分误差不可控。失败路径保证调用者 state 与 endpoint 的最后可信状态一致。 - directed inside→outside 事件的 DP 子积分定位使用同一 DP 误差控制（不混用 RK4）， 成本计入 endpoint 但不伪造 accepted 主步。定位的 time bracket 停止容差取局部步长量级与真实 `nextafter` ULP 的组合（不随远端 `|t|` 放大），任何 可表示的非零 target−before 区间都实际积分，只有 target==before 才允许直接取状态； `t+h==t` 明确失败或返回已接受可信状态，不伪造轨迹。子积分失败经独立 RayReason 透传 （`TIME_UNAVAILABLE→TIME_RANGE_EXHAUSTED`、`INVALID_METRIC`、`OUT_OF_DOMAIN`、 `INTEGRATION_ERROR`），与 worldtube descriptor 自身的 protocol/history 失败区分。 - 嵌入误差估计不保证找到全部几何事件；窄进出、first-crossing 与同一步内多个事件的 顺序必须独立验证，不能只凭步端点符号或误差通过认定事件完整。积分精度还须通过 逃逸方向、频移、null 残差及适用 backend 的独立守恒量验证；不靠反复归一化掩盖误差。 - DP54 accepted 区间的四次 dense interpolant 用于枚举事件候选。对分段匀速球形 worldtube，代入 `F=|x-center|²-radius²` 得到至多八次多项式；通过导数根划分 单调区间并隔离根，不能仅查看主步两端符号。根按过去传播方向排序，区分真正的 inside→outside crossing 与切触；motion 段边界截断主步并重新建立候选。 相机相对能量阈值的候选从同一 accepted 区间的 `L-L0` 获得。候选必须经误差受控 子积分验证，escape 与暗阈值按可信轨迹顺序处理；数值不可分辨的同时事件采用明确、 一致的优先规则，不能由 end descriptor 的排列决定终态。 `L-L0` 在动态时空中不假定单调：即使主步两端均低于阈值，也要处理步内的首次 upcrossing。阈值根按导数根分段，在局部 bracket 内确认并定位；不能用整个区间的 二分取代首次 crossing。所有候选的可信时间须在子积分后重新比较；当前不可分辨的 同时事件取 escape 优先。阈值是 `>=` 的闭事件，真正的步末 crossing 不能因根 合并到端点而丢失。无法确认事件时有界缩步，耗尽后报告积分失败，不回退为假成功。 `LOG_P0` 辅助监测策略保留 accepted-state 检查，不提供这一多项式事件定位能力； 生产 CLI 使用相机相对 `L-L0` 策略。 - lens-map 的积分 provenance 必须保存稳定的 stepper wire code、分组绝对/相对 容差、初始步长及上下限、拒步上限、初始计算/回溯配额与重试政策。render-only replay 消费文件策略，不能由当前 CLI 默认值覆盖历史配置，也不从缺失字段推断 自适应积分来源。未知 stepper code 或无效配置须拒绝；只含固定 RK4 来源的旧文件 可按明确的固定步语义导入，但不能补成 DP54 结果。 - accepted、rejected 与实际 RHS 成本按持久 sample 身份统计；续追保存累计值，回填 替换而不是再次累加。off-mesh witness 与原地提升后的顶点是同一样本，不能重复 计费；事件子积分和失败 stage 的 RHS 也计入真实成本。渲染地图保存这些诊断值， replay 的来源与成本统计应与原始 trace 一致。 --- # 19. 恒星 catalog 的内部表示 不要预存 RGB。 恒星应保留： - 无穷远天球方向； - 温度 \(T\)； - 辐射 normalization / amplitude。 概念结构： ```c typedef struct { float n[3]; float temperature; float amplitude; } Star; ``` 原因： ray tracing 会给每个 image 一个频率比： \[ g= \frac{\nu_\text{obs}} {\nu_\text{emit}} \] 对于黑体谱： \[ T_\text{obs}=gT_\text{emit} \] 因此温度是更自然的源数据。 最终根据： - intrinsic temperature； - frequency shift \(g\)； - lensing magnification \(\mu\)； - source amplitude； 得到 observed blackbody spectrum，再转换为 linear RGB。 --- # 20. 点源渲染 不将恒星 catalog 插值为天球亮度纹理。 最终 frame mesh 已知局部： \[ (x,y) \to \hat n_\infty \] 对每个局部可逆 triangle，反过来得到： \[ \hat n_\infty \to (x,y) \] 并查询 source triangle 内的 catalog stars。 ### 固定 2MASS tile 的基准结论（2026-08-27） 在 `5d80bb1` 上，以 54,160-source 的 2MASS `1 deg` tile、完全位于该 tile 内的 `0.6 deg` Minkowski 视场、1920 x 1920 输出和 **6 个 catalog render worker** 测得：每个 source triangle 扫描完整 tile 的有效墙钟成本约为 `105.9 us`，并且 wall time 对 triangle 数的线性拟合为 `11.855 s + 105.88 us * T`（`R^2 = 0.99955`）。完整记录、命令和原始时间见 [`benchmarks/2mass_fixed_tile_triangle_scan_2026-08-27.md`](benchmarks/2mass_fixed_tile_triangle_scan_2026-08-27.md)。 这说明局部稠密 tile 的全扫本身可以接受；catalog spatial index 的首要职责是 让每个 source triangle 只枚举相交的固定 1-degree tile，而不是扫描全天。不要 仅为把叶子降至数百颗星而预先激进细分；是否需要对极稠密 tile 再做浅层细分， 应由实际 Schwarzschild/BBH source-triangle 覆盖和负载 profile 决定。 每颗星获得： - image-plane 连续位置 `(x,y)`； - frequency shift \(g\)； - magnification \(\mu\)。 然后： \[ T_\text{obs}=gT \] 最终用一个 PSF： \[ P(x-x_s,y-y_s) \] 直接 splat 到 HDR framebuffer。 不要先把星压进最终 pixel 再卷积，因为那会丢失亚像素位置。 `--fast-mode` 是显式预览近似：把每个星像作为 delta 累积进 `N×N` 超采样 buffer（`--fast-deposit nearest` 只写 1 个超采样像素，`bilinear` 写 4 个 相邻像素），最后对整张 buffer 用单个全局离散核做一次卷积，再以 `N×N` box 平均降采样。核取目标 Moffat 在超采样尺度（宽度 `Nα`、同一 `beta`）下 的 pixel-area 积分且不做额外归一化；`1/N²` 的 box 平均正好重建最终 pixel-area 积分，所以核本身不改变指定的 FWHM/beta。`nearest` 在吸附到的 超采样格点上精确保持 PSF 形状，代价是至多 `1/(2N)` 个输出像素的位置量化； `bilinear` 精确保持亚像素质心，但会略微展宽 FWHM，因此默认是 `nearest`。 沉积误差测量、`N` 与取舍见 [`benchmarks/fast_mode_deposit_2026-09-18.md`](benchmarks/fast_mode_deposit_2026-09-18.md)。 该模式只支持 CPU PSF 后端；`--max-cache-psf-flux` 不适用（所有事件共用同一 个全局核半径）。CPU 构建中该全局卷积实现为零填充的 double 精度 FFTW 线性 卷积（`fftw3`/`fftw3_omp`）：核频谱、FFT plans 与 planar scratch 在初始化时 一次性建立并跨帧复用，`resolve()` 只做 zero/pack、批量 R2C、频域乘、批量 C2R，以及在 `(R,R)` 处的裁剪、`1/(Pwidth·Pheight)` 与 `1/N²` 归一化和 `+=` 到 HDR。旧的嵌套空间循环作为测试/基准参考保留，不是运行时 fallback。 依赖与 plan 模式取舍见 `build.md` 与 [`benchmarks/fast_mode_fftw_2026-09-25.md`](benchmarks/fast_mode_fftw_2026-09-25.md)。 plan 策略可用 `--fast-fftw-plan estimate|measure|wisdom|wisdom-update` 显式选择； `wisdom` 用 `FFTW_WISDOM_ONLY` 并要求 `FILE` 旁的同源 `.meta`（FFTW 版本、 double 精度、supersample、FFT 尺寸、核半径、worker 数）全部匹配，不匹配时明确 失败而不静默重 plan；`wisdom-update` 先 measure，再把 wisdom 与 meta 各自写入 唯一临时文件并分别原子 rename（meta 最后提交），二者配对不具事务性；任一写入、 导出或 rename 失败都会使初始化失败。 movie catalog 预取提升到 movie 级：所有 frame mesh 完成 refinement 后，把所有 可用 source triangle 涉及的 tile 合并为一个固定大小 `CatalogTileSet`，按批 （默认每批 `512` 个 tile）并行读取、每批末串行提交一次；逐帧 splat 只读已 immutable 的 cache，不再逐帧做 tile 扫描、prefetch OpenMP 区域或打印 prefetch 日志。单帧仍保留 frame 级预取；内存 catalog 完全不做 tile 扫描；imported multi-frame map 与 observer movie 走同一 union 路径。 movie PNG 编码/写盘由一个单 producer、单 writer 的有界队列（默认容量 2）承担， 与下一帧渲染重叠。producer 完成 sensor bloom、tone mapping 和 sRGB 转换， 生成 clean RGB8；可选诊断网格由去重边光栅化为预乘 alpha 的 sRGB RGBA8 层。 HDR 与临时线段由 producer 释放；job 独立拥有这两个输出 buffer，成功 enqueue 后 ownership 转交队列。writer 先写 clean 图，再以 source-over 原地合成诊断层、写出 mesh sibling，最后释放 job buffer。单帧、movie 与 lens-map replay 共用此合成规则。 writer 的首个错误持久保存，使后续 submit 立即失败；`finish()` drain 已接受 job 后 join writer，所有退出路径都必须 join，绝不为求重叠而提前打印 `Rendered ... ok`。 诊断网格是最终 mesh 的只读可视化：一像素抗锯齿边按端点终态分别着色相邻半边， 在中点切换颜色。合成位于 tone mapping 与 sRGB transfer 之后，clean 图与 HDR 保持独立。调色与透明度配置见 [`usage.md`](usage.md)。 fast mode 的单星精度由 deposit 模式与 `N` 决定：`nearest` 的格点间距是 每轴 `1/N` 个输出像素，单帧瞬时舍入误差至多是 `1/(2N)`；在 `--psf-fwhm-pixels` 不变时它与图像分辨率无关，只有增大 `N`（或使用保持 质心、但会展宽 profile 的 `bilinear`）才会降低。 但单星偏移不等于一张稠密静态图的整图误差。生产密度 4K 银心单帧对照实测 fast 对标准路径的 HDR relative L2 为 `N=2` 时 `10.44%`、`N=4` 时 `4.37%`， 线性 RGB 总和偏差 `+0.018%`。把 FITS 三通道展平后，按 `|fast-base|` 排序 的前 `0.01%` RGB channel samples 占平方误差能量的 `95.7%`（`N=2`）或 `94.5%`（`N=4`），其 baseline linear-RGB 总和约占 `9%`。不重叠的 base-value 分区中，承载 flux 的暗值区间（`[1e-3,0.06)`、`[0.06,0.1)`、`[0.1,0.23)`， 合计约 `47%` 的 linear-RGB 总和）给出约 `2-5%` local relative RMS 与不高于 `0.19%` 的 signed bias。这一分布与大量暗星 PSF 的位置误差部分相消、已分辨 亮星核心保留亚像素位移的解释一致。因此只能说这次高密度 4K 单帧的银河纹理 明显好于全局 L2 数字所暗示，不能据此外推时间连续性或逐像素精度。 对 movie，平滑移动的星像跨越 nearest 格点边界时会跳变完整的 `1/N` 输出 像素：`N=2` 每轴 `0.5 px`、`N=4` 每轴 `0.25 px`。即使在 4K，这种不连续 跳变也可能形成可见闪烁，因此 nearest fast mode 尚未获得保持视频视觉特征的 资格。`bilinear` 可消除这一质心跳变，但存在上述 profile 展宽，也尚未作为 最终 movie 路径验收。定位上仍是近似路径：对单帧稠密星场，`N=2` 是吞吐 默认、`N=4` 是更高保真预览；per-event 参考 PSF 路径仍是定量与 movie 基准。 单帧误差分解可用 [`scripts/fast_mode_error_decomposition.py`](scripts/fast_mode_error_decomposition.py) 复现。 PSF 第一版可用 Gaussian； 以后可换成 Airy 或其他相机模型。 --- # 21. Lens magnification 局部 mapping： \[ (x,y) \leftrightarrow \hat n_\infty \] 自然给出 Jacobian。 点源每阶像的总 flux 应乘： \[ \mu = \left| \frac{d\Omega_\text{image}} {d\Omega_\text{source}} \right| \] 可由： - triangle area ratio； - 或局部 interpolation Jacobian； 计算。 critical curve / caustic 附近需要更多 adaptive refinement。 --- # 22. 并行策略 ## 22.1 ray evolution 使用 OpenMP 对当前 active ray pool 做 coarse-grained parallel loop。 优先： ```c #pragma omp parallel for schedule(static, CHUNK) for (...) { evolve_ray_through_current_slab(...); } ``` 若 strong-field rays 工作量差异明显，再测试： ```c schedule(dynamic, CHUNK) ``` 当前 movie `RayPool` 使用 `schedule(dynamic, 32)`：samples 按 frame 排列，连续静态分块会把同一 time slab 的 active rays 集中到少数 worker。 32 条一块的动态分配让空闲 worker 继续领取工作，也允许不同光线的积分开销 不一致。每条 ray 仍由一个 worker 推进，pool slot 与 frame/sample ID 不移动； 没有新增共享可变缓存或 per-ray allocation。 这一 chunk 大小来自 [393 帧 Schwarzschild 对照实验](benchmarks/ray_pool_scheduling_2026-09-06.md)， 并非所有 backend 的最优参数。实验同时比较连续静态、轮转静态和动态分块， 用完整 lens-map 文件校验调度前后的结果一致性。当前仍遍历整个 pool； active-index compaction、activation 和 mesh analysis 的并行化属于后续优化。 不要： - 一 ray 一个 OpenMP task； - 每 ray malloc； - 细粒度 mutex。 --- ## 22.2 adaptive mesh analysis recursive refinement 逻辑用**迭代 work queue**实现，不用真正的递归 task。 概念： ```c while (!queue_empty()) { tri = queue_pop(); if (triangle_is_good(tri)) { accept(tri); } else { request_missing_samples(tri); } } ``` 并行粒度应是： ```text (frame, root_tile) ``` 或至少按 frame。 每线程使用自己的 local queue 和 local request buffer。 最后 merge： - sort； - deduplicate； - assign vertex IDs； - bulk append。 避免多个 threads 同时修改一个 global mesh。 --- ## 22.3 可选 HIP PSF 后端 CPU reference 验证后，HIP 只替换最终 PSF 累积；catalog/source-triangle 查询、 inverse mapping、频移和黑体积分仍由 OpenMP worker 按 triangle 粗粒度并行完成。 每 worker 独占最多 16,384 项的 event chunk，不创建 per-star task。GPU sink 独占 一份 immutable cache 和 double RGB HDR；仅在 chunk 提交和 direct fallback 边界 使用共享锁。direct fallback 的完成 GPU 工作、下载、CPU 累积、上传必须在同一锁内。 每个 PSF 用 32 个线程协作访问相邻像素，仍使用原来的四点 phase cache 插值、圆形 support 和 double atomic 累积。两个 pinned host/device staging slot 只有在先前 kernel 完成后才能复用；所有 GPU 工作保持单 stream 顺序，不要求整个帧事件驻留内存。 计时 event 随 slot 复用，覆盖所有 batch，不再只测前 64 批。 以上参数由 [RX 9070 有界实验](benchmarks/hip_psf_2026-09-07.md) 支持，尚不是跨设备 最优值。该实验未运行完整全天渲染，不能用缩小样本直接宣称整帧加速比。分桶规约、 降低 HDR 精度或近似黑体积分均未作为此次优化的一部分。 2026-09-07 后续实验提供了独立的有界事件捕获/回放与 producer 诊断工具，详见 [回放记录](benchmarks/hip_psf_replay_2026-09-07.md)。诊断钩子仅编入测试程序， 不进入 renderer；输入限制为小型 CSV、指定 triangle 范围与事件上限，按既有 catalog/inverse-map/频移/黑体/cache/direct/min-Y 路径准备事件，不能输出伪 HDR。 索引诊断将有限子集预载进现有一度 tile 数据结构，在 worker 启动前完成，保留只读查询。 实验回放与 production 32-thread/event kernel 并存，比较 16×16 与 32×32 像素归约。 CPU 按完整、经 cache 限制的 support 矩形建立 tile 引用；GPU 仍按原来的圆形行边界 过滤贡献。phase/行边界在每个 chunk 内预计算并复用；每个像素线程使用 double RGB。 单列表 tile 独占写回 HDR；超过 256 项的列表拆成部分和，随后按固定 task 顺序合并。 引用上限 32 MiB、task 上限 2 MiB、部分和上限 128 MiB；缓冲由回放 sink 独占， 单 stream 顺序复用，超限明确失败。wing-clipped 事件保留原物理 support，不能在读取时 强制截短；只有 cache 遍历范围受限。direct 边界完整执行 finish/download/CPU/upload。 上述参数是实验配置，不是新的生产默认值。是否集成必须看包含初始化、分桶、预计算、 合并、清理和下载的完整成本及分布依赖，不能只引用无 atomic kernel 的单项耗时。 本轮生产 HIP 路径保持原实现，尚未进行新的完整全天渲染。 2026-09-11 的后续有界实验加入了按 chunk 的分布选择器：先统计 32px 中心 tile 占用， 再在高中心密度 chunk 使用 16×16 像素归约，其余保留 atomic。RX 9070 的配对回放中， 它在银心密集输入接近固定 tile16，同时在同贡献分散输入选择 atomic，详情见上述回放记录。 此选择器仍仅存在于测试程序；阈值来自固定样本，尚未验证 production producer/sink 锁竞争 或完整全天事件序列，因此不能成为 renderer 默认策略。 2026-09-11 又加入了 `PSF_BACKEND=dummy` 的只统计后端。它要求导入 `.grlens`，复用 production 的 OpenMP `schedule(dynamic,1)`、每 worker 16,384-event chunk，以及完整的 catalog/inverse-map/颜色/flux/cache 分类，但不分配完整 HDR，也绝不写 PNG/FITS。对银心 45 度 Schwarzschild 全天输入的一次 16-worker 全量运行中，598,264,924 个 cache event 形成 36,508 个满 chunk；满 chunk 跨越的 event-producing triangle 数为 p50=2、p75=3、 p90=10、p99=31、max=59，32px 中心 tile 占用为 p50=1、p75=2、p90=5、p99=30、 max=45。全部满 chunk 都满足当前测试选择器的 tile16 条件。 这个结果支持“真实宽视场下一个 worker 的满 chunk 通常集中在少数三角形/少数局部 tile”这一调度假设，但也揭示了稀有的远距离簇：满 chunk 的三角形质心包围盒对角线 p99 约 3744px，最大相邻三角形跳跃 p99 也约 3744px。因而不能把“少数三角形”等同于 “全 chunk 单个连续矩形”；production tile 原型仍应按稀疏 occupied-tile 列表处理，并用 真实流式端到端数据验证分桶成本。完整命令、hash、输出保护和原始记录见 [dummy chunk 统计](benchmarks/dummy_psf_chunks_2026-09-11/README.md)。 2026-09-12 在 production HIP sink 中加入了显式可选的 adaptive tile16 路径。默认仍为 原来的 atomic；设置 `HIP_PSF_ACCUMULATION=adaptive` 后，每个至少 8,192-event 的 chunk 先统计 32px 中心 tile，占用密度达到每 tile 32 events 时才使用 tile16，否则继续 atomic。 tile 路径按完整 cache-limited support 构建稀疏 touched-tile 引用，让一个线程拥有一个像素， 在寄存器中完成 double RGB 累积后写回 HDR。超过 256 项的 tile 列表使用有界部分和； references、tasks、partials 的 device 上限分别为 32 MiB、2 MiB、128 MiB，任一容量不足时 该 chunk 明确回退 atomic。event、metadata 和 partial 缓冲由 sink 持有并复用；所有工作仍在 同一 stream 中保持顺序，direct fallback 仍使用完整 download/CPU/upload 边界。 RX 9070 上，固定 65,536-event 银心真实输入连续回放 128 轮时，production atomic/adaptive 完整 replay wall 为 7.164/4.631 s（1.55×），其中 GPU 阶段为 7.132/3.875 s；adaptive 另含 0.607 s CPU 分桶和 0.091 s 上传。同贡献分散输入会选择 atomic。三 tile renderer 小样中 adaptive 选择 14 个 tile chunk、2 个 atomic chunk，splat wall 为 0.236 s，对照 atomic 为 0.285 s；两张 PNG 的 SHA-256 相同。这些有界结果随后被完整全天测试否定：相同 598,264,924 events 和 36,530 batches 下，atomic/adaptive 的 splat wall 为 480.500/586.763 s，process wall 为 533.62/640.35 s。tile kernel 虽从 477.302 s 降至 416.265 s，但 support 分桶串行花费 155.302 s，且 36,515 batches 被中心密度规则选为 tile16。中心位置的 32px tile 占用不能代表 47px PSF support 扩张后的 references、tasks、 列表长度和成本；固定事件样本也没有代表全天 tile kernel 的工作量分布。atomic 继续作为 默认值；上述旧串行 adaptive 是历史负结果，不能作为新版并行预分桶的结论。当时提出度量完整 support 工作量，并先证明分桶可并行隐藏或移至 GPU，再考虑完整全天复测。原始日志和误差检查见 [HIP 回放记录](benchmarks/hip_psf_replay_2026-09-07.md)。 随后实现的 producer 私有并行预分桶原型保留 GPU HDR/stream 的单 sink 所有权，且 mixed direct fallback 与小型 renderer 的图像正确性通过。固定银心事件 512-batch 回放中，sink 内 串行分桶和 16 worker 预分桶的完整 wall 分别为 3.886/3.958 s，后者没有收益；此输入又明显 低估全天真实的分桶工作，因此当时证据不足以推荐 adaptive。新诊断记录完整 support 生成的 references、tasks 和 merges。下一步需先获取真实全天 chunk 的这些指标及并行阶段 wall，而不是继续调中心密度阈值或直接重复昂贵的全天 HIP 渲染。 `DUMMY_PSF_SUPPORT_STATS=1` 可在无 HDR/GPU 的相同 16-worker chunk 流上报告这些指标的 分位数及第一遍 support 计数耗时；该耗时不包含 references 写入、task 构建或上传。 完整全天诊断的 36,508 个满 chunk 中，16px support references 的 p50/p90 为 746,749/802,816，touched tiles 的 p50/p90 为 49/168，最大 tile 列表的 p50 已达整批 16,384 个事件；全部 chunk 第一遍计数 worker 时间之和为 29.347 s。旧固定银心回放的 references 总量相近，但展开到更多 tile、最大列表更短。即使构造恰好 49 tiles、最大列表 16,384 的反事实热点回放，tile GPU 阶段仍约 7.3 ms/batch，远低于真实全天均值 11.40 ms/batch；15 个忙等 CPU 线程也未复现该差距。真实全天进度显示中段连续数万批 约 12–13 ms/batch、首尾约 7–8 ms/batch，因此尚不能把性能损失归因于某一个聚合几何量 或单纯主机抢占。当时提出采集旧版中段真实 chunk；该追溯现已搁置，不作为新版推进前提。 2026-09-13 用户重新构建并行预分桶 adaptive 后完成全天渲染：splat 282.944 s、GPU event 阶段 256.184375 s、完整进程 335.75 s。stash 回 `2d01c2c` 重编译后的原 atomic 为 481.417/478.336515/534.37 s，与历史 atomic 相差约 0.2%；该 HEAD 忽略 adaptive 环境变量。 相同 598,264,924 events 下新版 splat 快约 1.70×，但新版全天 HDR 尚未在本次文档工作中 比较。新版 bin 134.374 s 是 worker 累加时间，不能作为串行成本；旧 adaptive GPU 阶段 416.265 s 到新版 256.184 s 的差异仍未定位。**保留历史证据，搁置旧 adaptive 专项追溯， 后续以新版为优化基础；atomic 默认和参考路径本次不变。** 下一步目标是相同物理/数值参数、输入与 HDR/PNG 输出下完整进程达到 240 s/张以内， 这是待验证目标而非承诺。用户在 btop 连续观察到新版 GPU 约 94%、反复降至约 80%， CPU 满载，并报告原 atomic 能吃满 GPU；旧串行 adaptive 未观测，不能推断其利用率。 GPU busy 不代表净算力利用率，CPU 忙等需通过锁等待和线程 CPU time 定量确认。 仅消除当前 splat 与 GPU event 区间的约 26.8 s 差距仍需约 309 s/张；若其他阶段约 52.8 s 不变，目标要求 splat 约 187.2 s，需要进一步降低 GPU 执行或区间内间隙成本。 先利用已有全天输出完成 HDR 正确性比较，再在新版有界真实流式样本上拆分分桶、锁等待、 GPU 完成等待、staging/上传与 prepare/accumulate/merge，辅以同步设备遥测。根据证据 评估双缓冲重叠准备/上传，必要时对照专用提交线程和有界队列；均为待实验方案，保留单 sink HDR ownership、完成后复用、有界背压及 direct fallback 顺序。随后按实测占比优化 主要 kernel 或输入阶段，保留 double RGB 和现有误差阈值。先完成有界配对、正确性及 完整 provenance，再另行安排全天确认，不以恢复旧版慢因或 btop 达到 100% 为验收条件。 具体目标、计时口径、用户原始证据与限制见 [2026-09-13 全天复测与下一步目标](benchmarks/hip_psf_replay_2026-09-07.md#2026-09-13-全天复测与下一步目标)。 2026-09-12，黑洞渲染的黑体颜色计算改为以仓库内数值 LUT 为默认实现，取代原先 `optics.c` 中的 Wyman 解析 CIE 拟合。正常区使用仓库唯一的 `assets/blackbody/cie1931_2deg_xyz_1024.grbblut`：以 40 位精度把 Planck 谱和 `assets/CIE_xyz_1931_2deg.csv` 的 360--830 nm、1 nm 分段线性 CIE 数据积分，生成覆盖 `[670.146556,101408.88] K` 的 1024 个 `log(T)` 等距 XYZ 节点。`GRBBLUT3` header 固定标识该 reference，loader 校验 magic/version/endianness/物理常数/范围/长度与 payload checksum，并拒绝旧的 Wyman/5 nm 表与 runtime 造表。表内使用四点三次 Lagrange 插值， 在全部 1023 个区间中点相对独立 Gauss-16 积分的最大 XYZ channel 相对误差为 `5.4633e-5`，没有负 XYZ。由 1 nm 线性插值自身约 `3.2e-4` 的最坏误差估计，正常区采用 `1e-4` LUT 内部限与 `5e-4` 整体误差预算，避免继续追逐缺乏物理意义的 `1e-6`。 两端使用同一 CSV reference 的表外模型。高温端使用连续积分得到的 `A T+B+C/T` Rayleigh--Jeans 展开，在 `[Tmax,10Tmax]` 最大 XYZ 相对误差为 `1.6635e-5`。 低温端针对 X/Y 的 830 nm 端点和 Z 的 650 nm 零端点分别展开：在 `4.097321 K` 以下使用 crossover 锚定、在 `T -> 0` 趋于极限的端点式，其上用六段八次 `1/T -> log(scaled XYZ)` Chebyshev--Lobatto 插值接到 LUT。低温端点段/过渡段的最大 XYZ 相对误差分别为 `6.69e-5`/`8.16e-5`，因而三段均满足 `1e-4` 内部限；log 输出避免把极端 低温的线性下溢误解为拟合硬下限。 最终 legacy Wyman/5 nm 与 CSV LUT 的 1920x1080 double-HDR 比较显示，每通道峰值归一化 最大绝对变化为 R `9.00e-3`、G `7.02e-3`、B `8.12e-4`；这是 reference 升级本身造成的 预期变化，而不是 LUT 插值误差。旧 Wyman/5 nm 积分实现、可切换的 integral/fixed/LUT A/B 后端、有界性能证据与观测温度采样器保留在实验分支 `codex/blackbody-cost-experiment`（tag `blackbody-cost-experiment-2026-09-12`）。 --- # 23. Metric evaluator 的 thread-local workspace nmesh 单点查询需要： ```text global x → locate AMR leaf → global/local coordinate inversion → Lagrange interpolation → derivatives ``` 应给每线程一个 workspace： ```c typedef struct { uint32_t last_node[...]; double lx[12]; double ly[12]; double lz[12]; /* scratch arrays */ } MetricWorkspace; ``` 连续 RK stage 中 ray 大概率仍在同一 node，因此优先： ```text is x still inside last_node? yes: reuse no: tree / neighbor lookup ``` 不要在 evaluator 内维护共享 mutable cache。 --- # 24. 顶层主要数据结构 ```c typedef struct ObserverTrack ObserverTrack; typedef struct SpacetimeSource SpacetimeSource; typedef struct MetricSlab MetricSlab; typedef struct RayPool RayPool; typedef struct StarCatalog StarCatalog; typedef struct { double t_camera; ObserverState observer; FrameMesh mesh; HDRImage image; } Frame; typedef struct { Frame *frames; size_t nframes; } Movie; ``` --- # 25. 顶层执行逻辑 概念 `main()`： ```c int main(void) { ObserverTrack observer; SpacetimeSource spacetime; StarCatalog stars; Movie movie; load_observer_track(&observer, ...); open_spacetime(&spacetime, ...); load_star_catalog(&stars, ...); init_movie(&movie, &observer, ...); movie_init_coarse_meshes(&movie); for (;;) { SampleRequest *req = NULL; size_t nreq = movie_collect_refinement_requests( &movie, &req ); if (nreq == 0) break; RayPool rays; rays_build_generation( &rays, &movie, req, nreq ); rays_trace_generation( &rays, &spacetime ); movie_install_ray_results( &movie, &rays ); free_ray_pool(&rays); free(req); } movie_render_star_catalog( &movie, &stars ); movie_write_output(&movie); return 0; } ``` --- # 26. `rays_trace_generation()` 的核心流程 ```c void rays_trace_generation( RayPool *rays, SpacetimeSource *st) { while (spacetime_has_previous_slab(st)) { MetricSlab *slab = spacetime_load_previous_slab(st); activate_rays_for_slab( rays, slab ); #pragma omp parallel { MetricWorkspace ws; metric_workspace_init(&ws); #pragma omp for schedule(static) for (size_t i = 0; i < rays->n; ++i) { if (!ray_is_active(rays, i)) continue; evolve_ray_to_slab_left_boundary( rays, i, slab, &ws ); } metric_workspace_free(&ws); } compact_ray_pool(rays); spacetime_free_slab(slab); } } ``` --- # 27. repo 的几个主要模块 这里讨论的是代码职责，不强调具体文件树。 ## `observer` 负责： - 读取 worldline； - 插值 observer event； - 插值/正交化 tetrad； - 给每个 frame 提供 observer state。 不负责： - geodesic； - metric； - image rendering。 --- ## `movie / frame` 负责： - frame 时间表； - 每帧 adaptive image mesh； - vertex / triangle； - refinement passes； - ray sample request； - 安装 ray endpoint results。 --- ## `ray` 负责： - 批量 ray state； - ray 初始化； - active/inactive 状态； - stream compaction； - endpoint 保存。 --- ## `geodesic` 负责： - 3+1 null geodesic RHS； - frequency shift evolution； - ODE stepper； - 在给定 metric slab 中推进 ray。 不负责： - disk I/O； - frame mesh； - star rendering。 --- ## `spacetime` 负责统一 backend： - Minkowski； - analytic Schwarzschild； - nmesh 4D numerical spacetime。 负责： - time slab； - metric interpolation； - spatial derivatives； - 必要的 temporal derivatives； - metric/data 状态（成功、时间不足、空间域不足、invalid metric、内部错误）； - 渐近端声明与 escape worldtube 能力。物理暗终态由能量阈值 policy 判定，不由 位置分类决定。 --- ## `nmesh_backend` 负责： - 读取 nmesh node 文件； - 从 filename 恢复 AMR tree； - node spatial lookup； - global/local coordinate conversion； - Lagrange interpolation； - Lagrange derivatives； - temporal interpolation； - puncture trajectory。 --- ## `catalog` 负责： - 读取 Gaia / 2MASS 等 catalog； - 转换为内部 `(direction, temperature, amplitude)`； - 天球 spatial index； - triangle/source-region 查询。 --- ## `optics / psf` 负责： - blackbody frequency shift； - temperature → linear RGB； - flux normalization； - PSF splatting； - HDR accumulation。 --- ## `output` 负责： - HDR frame； - exposure； - 线性 HDR → 可选 display operator → sRGB encoding； - tone mapping； - PNG/EXR； - 最终视频编码接口。 当前实现决策：预览用的 tone-mapped PNG/PPM 路径为 线性 RGB → display operator → sRGB 传输曲线 → 8-bit。默认 display operator 是逐通道 soft clip \(T_2(x)=\sqrt{\tanh(x^2)}\)（`--tone-map softclip --tone-map-p 2`）；历史 Reinhard \(x/(1+x)\) 以 `--tone-map reinhard` 保留，用于复现旧输出。FITS `--hdr-output` 绕过 tone mapping，始终写 pre-tone-map 线性 RGB。这里的 operator 只是当前简化显示管线，不是最终相机/ 传感器模型。 可选的后处理式 sensor bloom（`--sensor-bloom-limit E --sensor-bloom-transfer e`，默认关闭）在 PSF 累积完成的线性 HDR 上、display operator 之前运行，CPU 与 HIP、标准与 fast-mode、直接 trace 与 imported lens map 共用同一入口。每个 RGB 通道独立、各向同性地把超过有限响应上限 \(E\) 的 溢出按固定比例 \(e\) 传播到 8 邻域，其余 \((1-e)D\) 被 drain 吸收，指向图像外的 部分在边界丢失。它是唯象的有限响应近似，不模拟特定 CCD/CMOS/Bayer 结构； \(E\) 使用 exposure 之后的 renderer-scale 线性 HDR；\(e\) 同时决定每轮保留传播的 比例和有效传播距离；模型允许信号损失，不守恒。原始 `--hdr-output` FITS 在模型 运行前写出，因此始终是 bloom 前的 PSF HDR；tone-mapped PNG/PPM 与视频帧在模型 之后写出。视觉式多尺度 bloom 不在当前范围内。mesh 诊断层在 tone mapping 与 sRGB transfer 后合成。 --- # 28. 开发顺序 ## Phase 0 — 平直时空 + 星表 目的： - 验证 observer/camera projection； - 验证天球方向； - 验证 Gaia/2MASS catalog； - 验证 temperature/amplitude； - 验证 PSF； - 验证 adaptive image mesh 数据结构。 平直时空： \[ (x,y)\leftrightarrow \hat n_\infty \] 可以作为严格 benchmark。 --- ## Phase 1 — analytic Schwarzschild 建议使用 horizon-penetrating coordinates。 验证： - 红移暗阈值截断与 shadow； - Einstein ring； - multiple images； - adaptive refinement； - local inverse maps； - magnification； - frequency shift。 --- ## Phase 2 — numerical single Schwarzschild 使用当前已能跑的 FOCCZ4 moving-puncture Schwarzschild。 特点： - 开始存在 gauge relaxation； - 随后 metric 趋于稳定； - 可真实测试 4D metric streaming； - 可和 analytic Schwarzschild 做 regression comparison。 比较： \[ \Delta \hat n_\infty, \quad g, \quad \text{end/outcome classification} \] --- ## Phase 3 — FOZ4c BBH 等 FOZ4c： - 双曲性； - constraint damping； - 数值实现； - BBH 稳定演化； 成熟后接入。 renderer 顶层架构原则上不应为 BBH 重新设计。 --- # 29. 尚未决定的问题 需要以后通过 prototype / convergence test 决定： - 时间输出 cadence； - time interpolation order； - time slab 最佳内存大小； - BBH production node 数； - 4D 输出总数据量； - 相机相对能量阈值 `L-L0` 的默认值标定（当前 `L_dark=8`，可 CLI 覆盖）； - adaptive triangle refinement criterion； - critical curve 附近最大 refinement level； - Gaia 与 2MASS 的最终组合； - PSF 模型； - ODE integrator； - observer trajectory 标准； - 当前预览/PNG 默认采用逐通道 \(T_2(x)=\sqrt{\tanh(x^2)}\) soft clip， Reinhard 作为兼容模式保留；最终 production color management、传感器模型、 曝光标定和 HDR 视频编码规则仍未决定； - 是否需要 diffuse Milky Way background； - 阈值监测量当前统一采用相机相对增长 `L-L0`（对全部 backend 生效），不再使用 绝对 `L`、`ln(p^0)` 或 Killing 相对量；这些量均不得与真正的 infinity `g=E_camera/E_source` 混同。 --- # 30. 当前最核心的设计原则 整个 repo 可以压缩成以下几条： 1. **最终目标是动态 4D NR 视频渲染，不是单帧 ray tracer。** 2. **所有帧的 rays 按 coordinate time 整合，通过 time slabs 倒序推进。** 3. **metric backend 可替换，但顶层 scheduler 与具体时空无关。** 4. **相机是 worldline + tetrad track，不是固定三维坐标。** 5. **恒星是点源 catalog，不是 RGB 天球纹理。** 6. **恒星存温度与 amplitude，频移通过温度缩放自然处理。** 7. **image-plane 使用 adaptive local triangles 构造局部 inverse lens maps。** 8. **adaptive refinement 采用多 pass，避免和 out-of-core metric streaming 冲突。** 9. **ray evolution 用 SoA + OpenMP bulk loops。** 10. **recursive refinement 用 thread-local iterative queues，不用细粒度 OpenMP tasks。** 11. **nmesh 原生 DG node 输出直接作为空间插值数据，不做 Cartesian resampling。** 12. **时间导数优先由 3+1 identities、gauge RHS 或 temporal interpolant 解析获得。** 13. **所有高开销 mutable cache 都 thread-local。** 14. **第一版 CPU-only，先把物理与数据流做正确，再谈 GPU。** 15. **从 Minkowski → analytic Schwarzschild → numerical Schwarzschild → BBH 逐级验证。** ## Schwarzschild 自由落体轨迹辅助脚本 `scripts/schwarzschild_camera_track.py` 在入射 Cartesian Kerr–Schild 坐标中使用 $g_{\mu\nu}=\eta_{\mu\nu}+(2/r)\ell_\mu\ell_\nu$、$\ell_\mu=(1,x_i/r)$， 固定 $G=c=M=1$ 与解析 backend 一致。解析微分 metric 构造四维 Christoffel， 以本征时联合积分 $dz^\mu/d\tau=u^\mu$ 与 $de_{(a)}^\mu/d\tau=-\Gamma^\mu_{\alpha\beta}u^\alpha e_{(a)}^\beta$。 $e_{(0)}=u$ 同时满足自由落体方程；四加速度为零时费米–沃克输运就是平行输运。 初始坐标速度及指向构造沿用 10.1 的约定，也接受经验证的完整标架。 不在积分期间重新正交化以隐藏误差；输出正交归一误差超过 $10^{-6}$ 时失败。 DOP853 默认 rtol=$10^{-10}$、atol=$10^{-12}$，可配置并通过解析径向自由落体、 圆轨道和圆轨道平行输运的收敛回归验证。视界不终止相机；默认 $r=10^{-3}M$ 只是可配置的奇点数值保护边界，不等于精确撞击奇点。相机轨迹与光线终态解耦： 相机合法性与位置 cutoff 无关，光线正常暗终态由 §18 的能量阈值 policy 决定。 CSV 保持 21 列不变，记录 $\tau=k/\mathrm{fps}$ 与积分得到的真实坐标时间。 只保留不超过请求持续本征时或提前终止时刻的规则采样，包含 $\tau=0$。 `--movie-track-samples` 选择每个 CSV 样本直接生成一帧，绕过坐标时间均匀插值， 忽略 renderer 的 start-time/duration/fps。默认 movie 路径不变；两条路径随后共用 按真实 coordinate time 的 time-slab 调度，不把本征时冒充坐标时间。