逐像素弯曲时空：singularity-rs 黑洞渲染器深度解析

写下这个渲染器的过程，本质上是在和广义相对论的方程较劲——让它们在 GPU 上实时可解，又在物理保真与渲染预算之间走出一条平衡线。本文是对这次工程内核的解剖：不是使用说明，而是把每一条物理公式、每一处工程妥协都摊开来讲清楚。涉及到的数学，都尽量给出推导脉络，而不是只甩一个结论。

一、问题的本质：为什么黑洞渲染很难

渲染一个黑洞，难点不在于"画一个黑色的圆"——而在于光本身不再走直线。

在黑洞附近，时空被质量弯曲。一束从远方恒星射来的光，其轨迹会被引力折射，绕黑洞缠绕，甚至绕好几圈才到达你的眼睛。这意味着：屏幕上每一个像素的颜色，不能由一条直线射线与场景求交得到，而要由一条弯曲的测地线决定。这条测地线的形状取决于黑洞的质量、自旋，以及光线的撞击参数（impact parameter）——即光线若不受偏折时与黑洞中心的最近距离。

于是渲染退化为一个边值问题：对每个像素，从相机出发反向追踪一条测地线，在弯曲时空中数值积分，直到它（a）坠入事件视界——该像素为黑；（b）逃逸到无穷远——该像素采样恒星背景；或（c）途中穿过吸积盘——该像素叠加盘的辐射。

这个积分没有解析解（Kerr 时空的测地线方程不可积为初等函数），必须数值求解。而你要在 200 万像素上、每秒 60 次地求解它。这就是 singularity-rs 要攻克的根本矛盾：广义相对论的精度 vs. 实时光栅化的预算。

二、物理核心：从 Schwarzschild 到 Kerr 的精确退化

2.1 自然单位制：Rs = 1

整个项目贯穿一个约定：以史瓦西半径为单位。

史瓦西半径 R_s 是黑洞的"特征尺度"——对一个不旋转、不带电的黑洞，它就是事件视界的半径，定义为：

其中 $G$ 是引力常数，$M$ 是黑洞质量，$c$ 是光速。在数值实现里取 R_s = 1，于是由上式反推质量

这就是代码里到处出现的 let m = 0.5; 的来历。所有距离、速度、频率都以 R_s 为标尺。这不仅消除了量纲混乱，更让物理常数在 Rust 与 WGSL 两侧以同样的字面量出现——这是 CPU↔GPU 镜像一致性的基础。

在这个单位制下，史瓦西黑洞的临界撞击参数 b_{\text{crit}} 是一个关键量：它是刚好能"擦边"逃逸的光线的撞击参数，小于它则被捕获。精确值为：

这个值的来历值得展开。在史瓦西时空中，零测地线（光线）满足守恒的能量 $E$ 与角动量 $L$，撞击参数 b = L/E。结合测地线方程，光线的有效势能允许在 。结合测地线方程，光线的有效势能允许在 r = 1.5\,R_s 处存在一个不稳定圆轨道（光子球）。一条刚好在光子球外缘掠过的光线对应最大的可逃逸撞击参数，解轨道方程即得上式。

在渲染器里，b_{\text{crit}} 不是硬编码的——它从积分中自然涌现：撞击参数 不是硬编码的——它从积分中自然涌现：撞击参数 b < b_{\text{crit}} 的光线被黑洞捕获， 的光线被黑洞捕获，b > b_{\text{crit}} 的光线逃逸。但它被写成了 physics.rs 里的常量（因为 f32::sqrt 不是 const fn，无法在编译期算），并由测试守护：

pub const BCRIT: f32 = 2.598_076;

2.2 弯曲加速度：一条公式的两种命运

整个渲染器的物理灵魂，浓缩在 deriv 函数（着色器）/ kerr_bending_accel（CPU 镜像）里：

let radial = -1.5 * RS * h2 / r5 * pos;          // 径向弯曲（Schwarzschild）
let drag   =  2.0 * M * a / r3 * spin_axis.cross(dir);  // 参考系拖曳（Lense-Thirring）
let accel  = radial + drag;

先拆开看每一项。

角动量 $h$。 设光线位置 \vec{p}，方向 ，方向 \vec{d}（单位向量），定义

在球对称的史瓦西时空中，$h$ 就是每单位能量的轨道角动量，是一个守恒量。代码里 let h = pos.cross(dir); let h2 = h.dot(h); 即此。

径向弯曲项。 这一项把光线"拉向"黑洞中心，方向沿 -\vec{p}（指向原点）：

其中 r = |\vec{p}|。它的物理来历：光线的偏折可以由测地线方程的二阶近似得到，二阶导 。它的物理来历：光线的偏折可以由测地线方程的二阶近似得到，二阶导 \ddot{\vec{p}} 中主导项正比于 中主导项正比于 h^2/r^5 乘以位置向量。 乘以位置向量。r^5 = r \cdot r^4 在分母，保证了 在分母，保证了 r\to 0 时偏折急剧增强（强场）， 时偏折急剧增强（强场），r\to\infty 时迅速衰减（弱场近似牛顿）。代码里 r5 = max(r*r*r*r*r, 1e-6) 的 max 是防止除零的数值护栏。

参考系拖曳项（Lense-Thirring 前导项）。 当黑洞自旋 \chi = a/M \neq 0，旋转的质量会"拖动"周围的时空一起转。这一项沿自旋轴 ，旋转的质量会"拖动"周围的时空一起转。这一项沿自旋轴 \hat{a} 的横向方向施加偏转：

其中 a = \chi M 是 Kerr 自旋长度（几何化）， 是 Kerr 自旋长度（几何化），\hat{a} 是自旋轴方向（代码里固定为 是自旋轴方向（代码里固定为 +Y）。叉积 ）。叉积 \hat{a}\times\vec{d} 给出一个垂直于自旋轴和光线方向的横向分量，这正是"拖曳"的方向。 给出一个垂直于自旋轴和光线方向的横向分量，这正是"拖曳"的方向。r^3 在分母，意味着拖曳随距离立方衰减。

精确退化。 这个公式有一个极其优雅的性质：当 \chi = 0 时，a = 0，drag 项精确归零，积分器退化为纯粹的 Schwarzschild。这不是近似，不是 if 分支，而是数学上的精确退化。自旋参数 \chi 驱动着两个可观测效应：

1. 参考系拖曳（Lense-Thirring）：drag 项使光线的方向产生横向偏转，破坏了球对称，让光晕相对视线轴发生剪切。
2. 自旋依赖的捕获半径：见下一节。

两项的 \chi = 0 极限都精确回到 Schwarzschild——这被 14 个单元测试逐条断言。

2.3 自旋依赖的捕获半径：Kerr 视界

旋转黑洞的事件视界比史瓦西的小。Kerr 度规的视界半径 r_+ 为：

代入 M = 0.5、、a = \chi M = 0.5\chi：

• \chi = 0：：r_+ = M + M = 2M = R_s = 1（史瓦西）；
• \chi = 1（极端自旋）：（极端自旋）：r_+ = M + 0 = M = 0.5。

pub fn kerr_horizon(chi: f32) -> f32 {
    let m = 0.5;
    let a = chi * m;
    m + (m * m - a * a).max(0.0).sqrt()
}

.max(0.0) 是数值护栏：极端自旋下 M^2 - a^2 理论上为零，浮点误差可能让它略负。光线穿过 理论上为零，浮点误差可能让它略负。光线穿过 r < r_+(\chi) 即判定捕获。

2.4 ISCO：自旋如何塑造吸积盘的内缘

吸积盘的内边缘不是任意设定的——它追踪 Kerr 的最内稳定圆轨道（Innermost Stable Circular Orbit, ISCO）。在 ISCO 以内，物质无法维持稳定圆轨道，必然螺旋坠入黑洞，所以那里不再有盘的物质。

ISCO 没有简单的初等形式，但有 Bardeen-Press-Teukolsky (1972) 给出的闭式解。定义两个中间量：

则顺行（prograde，与黑洞同向旋转）轨道的 ISCO 为：

几个极限值便于校验：

• \chi = 0：：Z_1 = 3，，Z_2 = 3，代入得 ，代入得 r_{\text{ISCO}} = M(3+3-\sqrt{0}) … 化简后为 … 化简后为 6M = 3\,R_s；
• \chi = 1：：r_{\text{ISCO}} = M = 0.5\,R_s。

pub fn kerr_isco(chi: f32) -> f32 {
    let m = 0.5;
    let z1 = 1.0 + (1.0 - chi * chi).cbrt() * ((1.0 + chi).cbrt() + (1.0 - chi).cbrt());
    let z2 = (3.0 * chi * chi + z1 * z1).sqrt();
    m * (3.0 + z2 - ((3.0 - z1) * (3.0 + z1 + 2.0 * z2)).sqrt())
}

这里有一个精妙的工程设计决策：params.rs 里有一个 disk_inner 字段，但 mirror_params（每帧把参数镜像到 GPU uniform 的函数）故意覆盖它：

// disk_inner 是自旋导出的（Kerr ISCO）；params.disk_inner 字段被忽略。
u.disk_inner = crate::physics::kerr_isco(params.spin);

disk_inner 字段存在的唯一理由是 UI 默认读数的便利，运行时永远被 ISCO 公式覆盖。这消除了"用户手动设了个内缘但物理上不自洽"的整类 bug——单一数据源，物理保证。

三、积分器：GPU 上的自适应 RK45

3.1 为什么不用固定步长

黑洞附近的时空曲率变化剧烈：在事件视界附近，曲率极高，光线急剧偏转，需要极小的步长才能准确积分；而在远离黑洞的平直区域，大步长即可。固定步长面临两难——步长太小则在平直区浪费算力，步长太大则在强场区积分失准，捕获/逃逸判定出错。

答案是自适应步长：根据每步的局部截断误差估计，动态放大或缩小步长。

singularity-rs 选择了 Dormand-Prince RK45 方法——这恰好是 MATLAB ode45 背后的方法。它属于嵌入式 Runge-Kutta 家族，核心思想是用同一组中间求值同时算出两个不同阶的解，它们的差就是误差的"免费"估计。

3.2 RK45 的数学原理

考虑常微分方程初值问题 \dot{y} = f(t, y)。一个 $s$ 级显式 Runge-Kutta 方法形如：

其中各级 k_i 为：

系数 c_i、、a_{ij}、、b_i 构成 Butcher 表。Dormand-Prince 的精妙之处在于：它设计了两套权重 b_i（5 阶）和 （5 阶）和 \hat{b}_i（4 阶），共用完全相同的 6 个 k_i 求值，于是：

误差估计为两者之差：

这相当于"用 6 次函数求值的代价，同时得到了 5 阶解和 4 阶误差监控"——没有一次额外的求值浪费。

代码里就是这张 Butcher 表的逐字翻译。以第 3 级为例：

let p3 = pos + (k1p * 0.075 + k2p * 0.225) * dt;   // c3 = a31 + a32 = 0.3

其中 a_{31} = 0.075、、a_{32} = 0.225 是 DOPRI 表的精确系数。6 级、5 阶权重 是 DOPRI 表的精确系数。6 级、5 阶权重 (37/378,\,0,\,250/621,\,125/594,\,512/1771,\,0)、4 阶权重 、4 阶权重 (2825/27648,\,0,\,18575/48384,\,13525/55296,\,277/14336,\,1/4)——这些分数在 Rust 与 WGSL 两侧逐字一致。

3.3 步长控制律

有了误差估计 \text{err}，如何调整步长？标准做法是：希望误差接近（但不超过）容差 ，如何调整步长？标准做法是：希望误差接近（但不超过）容差 \text{tol}。对于 4 阶误差估计，局部误差按 。对于 4 阶误差估计，局部误差按 h^5 缩放（ 缩放（\text{err} \propto h^{p+1}，，p=4），所以最优步长满足：

令其等于 1，解出：

指数 1/5 = 0.2 正是代码里 powf(0.2) 的来历。实际实现加了上下限钳制和 $0.2$ 的指数：

dt = (dt * (tol / err.max(1e-12)).powf(0.2)).clamp(dt_min, dt_max);

这保证步长在平直区膨胀（上限 dt_max = 4×dt_init）、在强场区收缩（下限 dt_min = 0.25×dt_init）。

3.4 GPU 自适应的陷阱与破解

自适应积分在 CPU 上平平无奇，但在 GPU 上却暗藏杀机。GPU 着色器没有递归，循环边界有严格限制，而且绝不能让任何像素的死循环拖垮整个 warp。

march 函数的自适应循环有三个精心设计的语义：

语义一：预算只消耗"接受"的步。

var budget = steps_max;     // 硬上限（桌面 300）
loop {
    if (budget == 0u) { break; }
    let step = rk45_step(pos, d, dt);
    let err = step.err;
    if (err > tol * 10.0 && dt > dt_min) {
        dt = dt_min; continue;   // 拒绝：缩小 dt 重试（不消耗预算）
    }
    budget = budget - 1u;        // 接受：消耗一个预算单位
    ...
}

被拒绝的步（误差过大）缩小 dt 后重试，不消耗预算。这保证了"步数"参数真实反映了几何分辨率，而不是"尝试次数"。

语义二：dt_min 是强制接受的底板。

如果一条光线进入了某区域，那里的误差即使在最短步长下也超容差（极端自旋、视界附近），纯自适应循环会永远 continue 下去——GPU 挂起。dt > dt_min 这个守卫确保：一旦 dt 已经触底，无论误差多大都接受这一步：

if (err > tol * 10.0 && dt > dt_min) {  // ← 只有还能缩时才重试
    dt = dt_min; continue;
}
budget = budget - 1u;                    // 否则强制接受

这是一个防 GPU 挂起的安全阀。

语义三：步长按经典公式细化。 见 3.3 节。

这三个语义在 CPU 镜像 is_captured_rk45 里被逐字复制，并由 tests/physics_test.rs 的循环级测试守护——测试不仅验证单步，更验证预算、dt_min 底板、r_+(\chi) 捕获半径这些循环级不变量。

3.5 一个诚实的命名：RK45 是表，不是阶

代码注释里有一段罕见的诚实声明：

积分器在每个阶段后都对方向重新归一化（normalize）。这保持方向为单位向量（物理上，光线的切向量始终单位长），但破坏了纯 Dormand-Prince 表的严格误差阶。实测的误差衰减介于 2~4 阶之间，取决于几何。

作者没有粉饰这一点。自适应循环只依赖误差关于 dt 单调这一性质（rk45_step_error_shrinks_monotonically_with_dt 测试验证了它），而单调性成立。"RK45"命名遵循着色器的 Butcher 表，而非精度保证。这种对自身局限的坦诚，是工程成熟度的标志。

四、架构的支点：CPU↔GPU 物理镜像

整个项目最独到的架构决策，是将物理积分器在两个地方维护两份等价实现：

为什么要做镜像？ 因为 GPU 着色器无法单元测试。整个渲染器的物理正确性——捕获与逃逸的边界、自旋退化、ISCO 收缩——全部系于一段无法被 cargo test 触达的 WGSL 代码。如果物理只在着色器里，永远无法在 CI 里回归验证它。

镜像的设计哲学是：让物理可测的部分尽可能贴近不可测的部分。physics.rs 不是为了运行时性能（二进制几乎只用它的 kerr_isco/kerr_horizon 做 UI 读数），而是为了让捕获/逃逸边界这个最关键的物理不变量变得可测试。

代价是双份维护：改变一处物理，必须同步另一处。AGENTS.md 明确列出了必须锁步的函数清单（bending_accel / kerr_bending_accel / rk45_step / is_captured / is_captured_rk45）和三个循环级不变量（预算=仅接受步、dt_min 强制接受底板、r_+(\chi) 捕获半径）。Butcher 表的系数在 Rust 与 WGSL 两侧逐字一致。

这构成了一种人肉实现的同伦检查：两边做的事在语义上相同，测试验证 CPU 侧的语义，注释和纪律保证 GPU 侧跟随。脆弱吗？是。但目前没有更好的办法让 GPU 物理可回归。

五、渲染管线：七台相机的 HDR 合唱

物理积分只是第一幕。一个"物理正确但视觉平庸"的黑洞——纯黑圆盘加笔直光线——并不动人。singularity-rs 的第二幕是一条多阶段 HDR 后处理管线，把物理输出锻造成电影级画面。

5.1 全屏 quad 链

整个画面由七个 Camera2d 协作产出，按 Camera.order 从 −20（最先）到 composite（最后）严格排序：

offscreen quad ──→ Rgba16Float Image (HDR 场景)
    │
    ▼ order -19
brightpass  ──→ 半分辨率 HDR (亮度提取)
    │
    ▼ order -18 → -17
blur 下采样 ×2  (1/4 → 1/8 分辨率)
    │
    ▼ order -16 → -15
blur 上采样 ×2  (1/8 → 1/4 → 1/2 分辨率)
    │
    ▼ composite
composite   ──→ 窗口 LDR surface (ACES 色调映射)

每个阶段都是一个全屏 quad，挂载一个自定义 Material2d，读写 Rgba16Float（每通道 16 位浮点）的离屏纹理。

为什么是 16 位浮点？ 这是 HDR 的关键。普通的 8-bit sRGB 纹理每个通道只能存 $[0,1]$。但吸积盘内缘的辐射亮度可以远超 1.0（物理上，热等离子体的辐射强度没有上限），只有浮点目标能保留这个动态范围，供后续 bloom 提取和色调映射使用。Rgba16Float 的半精度（6 位指数 + 10 位尾数）足以覆盖渲染所需的亮度量级，又比 32 位浮点省一半显存。

5.2 Bloom：金字塔模糊

Bloom（辉光）模拟真实镜头/相机传感器对强光的溢出响应——亮的东西会"晕开"一圈光。管线采用经典的降采样-上采样金字塔。

Bright-pass 亮度提取。 从 HDR 场景中提取亮度高于阈值的像素。这里用了软膝（soft knee）而非硬截断：

let lum = dot(hdr, vec3(0.2126, 0.7152, 0.0722));   // Rec.709 亮度
let soft = max(lum - u.threshold, 0.0) / (lum + 0.0001);
let contribution = hdr * soft;

Rec.709 亮度公式 Y = 0.2126R + 0.7152G + 0.0722B 反映人眼对绿光更敏感的特性。软膝让近阈值的像素平滑过渡到零（而非硬切），避免提取出硬边。

金字塔模糊。 两次降采样（半→四分之一→八分之一分辨率）+ 两次上采样（八分之一→四分之一→二分之一）。每次卷积用 13 抽头的加权高斯近似核。为什么用金字塔而不是单次大模糊？因为低分辨率等价于大空间范围的模糊，但成本极低——在 1/8 分辨率上做 13 抽头卷积，等效于全分辨率上做约 100 抽头的大核卷积，但只花 1/64 的像素数。

上采样阶段以 blend 系数（0.6、0.8）与上一级混合，重建宽频率范围的辉光。最终 bloom_final 纹理在 composite 阶段叠加到场景上。

一个务实的简化：BloomQuality 枚举有 Off/Low/Medium/High 四档，但 spawn_bloom_pipeline 实际上总是发射完整的三级金字塔（High）。调用方仅以"是否调用"来区分 Off 与 On——部分金字塔留作未来增强。这是"先做对，再做好"的典型工程克制。

5.3 ACES 色调映射：从 HDR 到 LDR 的艺术

HDR 场景（亮度可达数十、数百）必须压缩到显示器能表达的 $[0,1]$ LDR 范围。直接线性缩放会让暗部全黑、亮部全白，丢失细节。色调映射曲线解决"如何优雅地压缩动态范围"。

singularity-rs 用 ACES（Academy Color Encoding System）的 Narkowicz 解析拟合：

fn aces_tonemap(x: vec3<f32>) -> vec3<f32> {
    let a = 2.51; let b = 0.03; let c = 2.43; let d = 0.59; let e = 0.14;
    return clamp((x*(a*x+b)) / (x*(c*x+d)+e), vec3(0.0), vec3(1.0));
}

这条曲线的形状值得品味：

• x \to 0（暗部）：（暗部）：f(x) \approx bx/e = (0.03/0.14)x \approx 0.21x——近似线性，但斜率小于 1，抬高了暗部，让暗处细节可见；
• x \to \infty（亮部）：（亮部）：f(x) \to a/c = 2.51/2.43 \approx 1.03——渐近于 1，形成柔和的胶片肩部（filmic shoulder），高光柔和滚落而非硬削顶；
• 中间段是平滑的 S 形过渡。

仅 5 次运算，却赋予了画面电影感的"胶片质感"。这是整个视觉"高级感"的来源之一。Composite 公式为：

六、吸积盘：从平板到体积的进化

吸积盘是黑洞渲染的视觉主角。它的难点在于：盘是一个三维的、有内部结构的、旋转的等离子体，而光线在盘内穿行时步长是自适应变化的。

6.1 体积积分与"辐条"bug

早期的实现用每步采样的密度乘以 RK45 步长来积分盘的不透明度：

这产生了一个诡异的视觉 bug：盘上出现径向辐条。

根因精妙。RK45 的步长是空间自适应的——黑洞附近步长小（dt_min 触底），远处步长大（dt_max），且相邻像素的步长网格不同。用 density × step_len 加权，等于把一个每像素不同的亮度调制沿径向注入——步长最小处（内盘）辐条最密。本质上，积分权重耦合了不该耦合的数值参数。

修复方案 integrate_disk_segment 展现了出色的工程洞察：将每步解析地裁剪到盘的几何板层（slab），按板层内的世界空间交线长度加权，而非按 RK45 步长。

数学上，设盘是 |y| \le H(r) 的板层（$H$ 是随半径变化的高度）。光线的这一步从 的板层（$H$ 是随半径变化的高度）。光线的这一步从 \vec{p}_{\text{prev}} 到 到 \vec{p}_{\text{new}}，参数化为 ，参数化为 \vec{p}(t) = \vec{p}_{\text{prev}} + t(\vec{p}_{\text{new}} - \vec{p}_{\text{prev}})，，t\in[0,1]。解 。解 |y(t)| = H 得到进入/离开板层的参数 得到进入/离开板层的参数 t_0, t_1，于是光线在盘内的几何长度为：

let seg_len = (t1 - t0) * length(new_pos - prev);
for (var i = 0u; i < N; i++) {
    let s = disk_color_volumetric(p, dir);
    let ds = s.density * seg_len / f32(N);  // 权重 ∝ 几何长度，非 RK45 dt
    ...
}

seg_len 只取决于光线入射角和当地板层高度，与积分步长完全解耦。辐条消失。

6.2 随机抖动消灭摩尔纹

即使修了辐条，盘上仍出现同心摩尔环。N 个采样点落在确定性的 (i+0.5)/N 网格上，与二维像素网格、与每条光线的 RK45 步格发生混叠——这是经典的采样混叠。

修复是分层随机抖动（stratified jittered sampling）。在第 $i$ 个分层（stratum，区间 [i/N, (i+1)/N]）内，采样位置加上一个该分层内的均匀随机偏移：

let stratum = f32(i) + hash13(vec3<f32>(pixel_seed, f32(i)));
let t = t0 + (t1 - t0) * stratum / f32(N);

关键性质：抖动幅度限制在分层内（\le 1/N），所以积分阶不变（仍是 $N$ 点求积），期望的积分值无偏。但确定性网格的相干摩尔带变成了非相干噪声，而 HDR + bloom 管线对非相干噪声的容忍度远高于相干带。

种子只用像素坐标（不含时间），于是噪声是空间稳定的——没有逐帧闪烁。这是两个层次的抗混叠智慧：物理积分层（解耦步长）和采样层（随机化网格），各自针对一种混叠机制。

6.3 黑体色与 Novikov-Thorne 温度

盘的颜色有两种模式。梯度模式是手工调的白热→深橙渐变加牛顿多普勒。黑体模式则走物理路线，由三部分组成。

(1) Novikov-Thorne 径向温度剖面。 标准薄吸积盘（Shakura-Sunyaev / Novikov-Thorne 模型）的局部有效温度随半径的分布为：

物理直觉：内缘（r \to r_{\text{ISCO}}）物质势能释放最剧烈，温度最高；外缘温度衰减。代码用 radial_temp 实现这一剖面：

let isco_r = clamp(inner / r, 0.0, 1.0);
let nt_factor = max(0.0, 1.0 - sqrt(isco_r));
let radial_temp = pow(isco_r, 0.75) * pow(nt_factor, 0.25);

(2) Tanner-Helland 黑体色温映射。 给定温度 $T$（开尔文），黑体辐射的峰值波长由维恩位移定律 \lambda_{\max} = b/T 决定（ 决定（b \approx 2.898\times10^{-3}\,\text{m·K}）。高温→短波→偏蓝；低温→长波→偏红。Tanner-Helland 给出了一个分段拟合，直接从 $T$ 算出 sRGB：

// 简化：T≤6600K 时 R=255，G/B 用对数拟合；T>6600K 时 B=255，R/G 用幂律拟合
let t = max(temp, 1.0) / 100.0;   // 归一化

6500 K 基准温度让内盘呈暖白色（接近 D65 标准白光），向外渐变为深橙——这正是 Gargantua 的标志性外观。注意输出要再转线性（pow(srgb, 2.2)）才能进 HDR 管线。

(3) Kerr 四速度多普勒因子。 这是黑体模式相对梯度模式最物理的地方。盘物质沿赤道圆轨道运动，其四速度 u^\mu 由 Kerr 度规的赤道分量决定。归一化条件 由 Kerr 度规的赤道分量决定。归一化条件 g_{\mu\nu}u^\mu u^\nu = -1 给出时间分量 给出时间分量 u^t，光子的能量比（观测者/发射者）为：

其中 k^\nu 是光子四波矢。代码里 kerr_doppler 实现了完整的赤道 Kerr 度规分量 g_{tt}, g_{t\phi}, g_{\phi\phi} 求解：

let g_tt = -(1.0 - 2.0*m/r);
let g_tphi = -2.0*m*a/r;
let g_phiphi = r*r + a*a + 2.0*m*a*a/r;
let u_t_sq = -(g_tt + 2.0*omega*g_tphi + omega*omega*g_phiphi);
let u_t = 1.0 / sqrt(max(1e-6, u_t_sq));
let l_photon = pos.z*dir.x - pos.x*dir.z;   // 光子角动量
let delta = 1.0 / max(0.01, u_t*(1.0 - omega*l_photon));

其中 \Omega = d\phi/dt 是开普勒角速度（见第七节）。这样得到的 是开普勒角速度（见第七节）。这样得到的 \delta 乘到温度上（ 乘到温度上（\delta 决定观测到的温度偏移），再用 决定观测到的温度偏移），再用 \delta^{3.5} 作 beaming（见 6.4）。接近侧蓝移变亮变蓝，远离侧红移变暗变红——这是梯度模式无法表达的物理特征。

6.4 相对论 beaming 的通用形式

盘的亮度和喷流的亮度都用到相对论 beaming。一个以速度 \beta = v/c 运动、方向 运动、方向 \hat{n}（指向观测者）的各向同性辐射源，观测到的辐射功率增强因子（多普勒增亮）为：

其中洛伦兹因子 \gamma = 1/\sqrt{1-\beta^2}。辐射强度按 。辐射强度按 \delta^p 增亮，幂指数 $p$ 取决于辐射机制：连续谱 增亮，幂指数 $p$ 取决于辐射机制：连续谱 p=3，某种几何平均取 ，某种几何平均取 p=3.5。

代码里盘用 \delta^{3.5}，喷流用 ，喷流用 \delta^{3.5}（并钳制上限 8×）。喷流速度 （并钳制上限 8×）。喷流速度 \beta = 0.92（接近光速的外流），这意味着接近的喷流那侧极其明亮、远离的那侧极暗——双极不对称。

6.5 相对论喷流：自旋驱动的极轴喷流

沿自旋轴（+Y）的双极喷流，关键的自洽性门控：

// 相对论喷流是自旋驱动的（Blandford-Znajek）：机制汲取能层，
// 而能层只在旋转黑洞存在。χ≈0 时没有驱动力。
if (uniforms.spin < 0.05) { return; }

Blandford-Znajek (1977) 机制表明，喷流的功率正比于黑洞自旋的平方 P_{\text{jet}} \propto \chi^2。即便 UI 开关 jets_enabled 为真，自旋为零时喷流也被物理地抑制——开关表达用户意图，自旋表达物理。没有这个门控，默认场景（spin=0）会在极轴上方出现毫无物理成因的蓝白光柱。

喷流的形态用高斯径向衰减 \exp(-r^2/\sigma^2) 和指数长度衰减 和指数长度衰减 \exp(-y \cdot 0.05)，加上沿轴向流动的湍流噪声（flow = pos.y*2 - time*8，让上下喷流分别向外流）。

七、轨道力学：Kerr 圆轨与 Lense-Thirring 进动

行星（如果开启）沿 Kerr 赤道圆轨道运行，并受 Lense-Thirring 效应发生轨道面进动。这给场景注入了动态的生命力。

7.1 开普勒角速度

赤道顺行圆轨道的角速度，Bardeen (1972) 给出：

pub fn kerr_orbital_frequency(r: f32, chi: f32) -> f32 {
    let a = chi * 0.5;
    1.0 / (r.powf(1.5) + a)   // M=0.5 时 √M=1/√2，此处做了单位归并
}

\chi = 0 时 时 a=0，，\Omega_\phi = M^{1/2}/r^{3/2} = \sqrt{M}/r^{3/2}，正是牛顿开普勒第三定律的形式 ，正是牛顿开普勒第三定律的形式 \Omega \propto r^{-3/2}——又一个精确退化。$a>0$（顺行）时分母变大，角速度比牛顿情形更小：自旋"帮"物质转，所以同样的半径转得更从容。

7.2 Lense-Thirring 节点进动

旋转黑洞拖曳时空，使得偏离赤道的轨道面会绕自旋轴缓慢进动——这就是 Lense-Thirring 效应。进动率 \Omega_{\text{LT}} = \Omega_\phi - \Omega_\theta，其中 ，其中 \Omega_\theta 是垂直 epicyclic 频率（Okazaki 1987; Kato-Fukue-Mineshige）：

let ratio = (1.0 - 4.0*a*m.sqrt()/r.powf(1.5) + 3.0*a*a/(r*r)).max(0.0);
let omega_theta = omega_phi * ratio.sqrt();
omega_phi - omega_theta

注释里特别强调了交叉项是 a\sqrt{M}/r^{3/2}（半径 $-1.5$ 次幂），不是 a\cdot\Omega_\phi/r——后者会让 ——后者会让 \Omega_\theta 偏大、进动偏小，且破坏"进动随 偏大、进动偏小，且破坏"进动随 \chi 单调增"的物理性质。这是一个容易写错、且错了不报错但物理不自洽的地方，注释的存在就是为了防止未来改坏。

\chi = 0 时括号 时括号 =1，，\Omega_\theta = \Omega_\phi，进动为零——史瓦西球对称，无进动。测试 nodal_precession_strong_field_exceeds_weak_field 还验证了强场区（$r<6$）精确进动率超过弱场近似 2aM/r^3：弱场近似只保留主阶，强场下高阶项变得重要。

由于真实进动极慢，UI 提供 planet_time_scale（默认 50×）放大模拟时间，否则进动在视觉时间尺度上不可见。

八、Bevy 0.19 的暗礁：五个会致静默失败的陷阱

这个项目最"接地气"的工程价值，在于它记录并解决了 Bevy 0.19 上一系列静默失败——不报错、不崩溃，只是画面空白或灰屏。这些 trap 的诊断过程本身就是宝贵的工程知识。

8.1 nudge_camera：静止相机不渲染

Bevy 0.19 issue #24448：静止的 Camera2d 在首帧后停止渲染。离屏相机是全场最静止的实体，它一旦冻结，合成相机会每帧重采样一张陈旧纹理——画面冻住。

修复极简而精妙：每帧给相机一个亚像素级的正弦抖动：

let nudge = (time.elapsed_secs() * 5.0).sin() * 1e-3;
t.translation.x = nudge;

10^{-3} 的平移远低于一个像素（世界单位下），视觉不可见，却让视图矩阵每帧变化，render graph 持续产出帧。振幅 的平移远低于一个像素（世界单位下），视觉不可见，却让视图矩阵每帧变化，render graph 持续产出帧。振幅 A\sin(\omega t) 的选择： 的选择：\omega=5 让运动平滑， 让运动平滑，A=10^{-3} 确保亚像素。这个 workaround 被明确标注"上游修复后移除"。

8.2 egui 上下文必须钉在合成相机上

bevy_egui 0.41 会把 PrimaryEguiContext 自动赋给第一个生成的相机——也就是离屏相机。但离屏相机渲染到 Rgba16Float 纹理，而 egui 的管线是 Rgba8UnormSrgb——格式不匹配崩溃。

修复链：disable_egui_auto_context（在 PreStartup 关掉自动赋值）+ 显式把 PrimaryEguiContext 组件加到合成（窗口）相机上。现在 bloom 管线有 7 个相机，这个陷阱尤其致命。

8.3 行星存储缓冲必须是真实的 ShaderBuffer 资产

Handle::default()（空句柄）会让 AsBindGroup 每帧返回 RetryNextUpdate，静默跳过 quad 的绘制——画面只剩相机的清除色（灰）。修复：启动时预填一个 MAX_PLANETS（32）大小的零缓冲，upload_planets 原地修改它，句柄永不变。

更微妙的是：绝不能每帧分配新的 ShaderBuffer。新资产还没有 GPU 资产，绑定解析又返回 RetryNextUpdate。必须原地修改已存在的资产，让 GpuShaderBuffer::prepare_asset 复用同一个 GPU 缓冲。

8.4 天空盒纹理绑定必须声明 dimension = "cube"

AsBindGroup 派生默认是 2D 纹理，但着色器声明的是 texture_cube<f32>。布局不匹配使管线特化失败，quad 静默不绘制。

8.5 WebGPU 上必须用 textureSampleLevel

天体盒在 RK45 主循环内被采样，而该循环控制流不均匀（if (accum_alpha > 0.99) { break; }，每像素提前退出）。WGSL 规范禁止 textureSample 出现在非均匀控制流（它需要屏幕空间导数做 mipmap 选择）。Chrome 的 Tint 强制执行并拒绝编译，而桌面端 naga 不执行——只在 Web 构建崩溃。textureSampleLevel 接受显式 LOD，允许在非均匀流中使用。

九、设计哲学的几条主线

通读全文，可以提炼出几条贯穿项目的设计原则。

9.1 单一物理数据源

disk_inner 被 ISCO 覆盖；捕获半径用 r_+(\chi) 而非硬编码；喷流被自旋门控。凡物理能决定的，不让用户/代码任意设。这消除了"看起来对但物理不自洽"的整类状态。

9.2 锁步四处的镜像纪律

镜像一个新参数到端到端，必须同步四处：BlackHoleParams（params.rs）→ BlackHoleUniforms（material.rs，注意 WGSL vec3 对齐填充）→ mirror_params 赋值（plugin.rs）→ 着色器结构体及使用。质量枚举还要加 UI 入口。这条纪律是物理正确性的程序性保障。

9.3 质量分级：桌面与 Web 的双轨

每个质量参数通过 cfg!(target_arch = "wasm32") 给出两套默认：steps 200/300、render_scale 0.5/0.75、bloom Low/High、disk Low/High、AA Off/Low。这不是简单的"Web 弱一点"，而是对每种质量维度独立权衡——RK45 积分已经是固定步长 RK4 的约 10 倍开销，所以 Web 的 render_scale 直接砍到 0.5 换取交互性。

9.4 调试设施内建，但克制

tests/physics_test.rs 守护物理不变量；preset_test.rs 守护参数预设；着色器内联所有函数（因为 naga_oil 跨模块导入曾静默失效）；注释详尽记录每个 workaround 的因果与上游 issue 编号。这些都是"为六个月后的自己/下一位维护者"写的。

9.5 资产嵌入：发布即自洽

Release 构建通过 bevy_embedded_assets（PluginMode::ReplaceAndFallback）把整个 assets/ 树嵌入二进制——发布的 singularity-rs 可执行文件是完全自洽的，旁边不需要 assets/ 文件夹。Debug 构建故意跳过嵌入：编辑 .wgsl 后 cargo run 立即反映，无需重编译 Rust。这个 debug/release 的区分，是渲染器开发循环的效率关键。

十、未竟之路：Phase 4 与精确 Kerr

当前实现用的是一个近似的 Kerr 弯曲加速度（径向 Schwarzschild 项 + Lense-Thirring 前导拖曳项），而非 Kerr 时空的精确 Cartesian 伪哈密顿量。精确形式涉及 Kerr-Schild 坐标下的伪哈密顿量，用 \Sigma、、\Delta 等度量函数表示：

精确方法还用 Carter (1968) 发现的第四守恒量实现 Hamilton-Jacobi 可分离性，能得到亚百分比精度的光子轨道。但在高自旋的 photon sphere 附近，近似方法的精度会下降——自旋依赖的捕获半径 r_+(\chi) 是精确的，但光线轨迹的偏折细节是近似的。

phase4-exact-kerr-experimental 分支尝试过精确 Kerr-Schild 哈密顿量：CPU 侧达到了亚百分比的 b_{\text{crit}} 精度（ 精度（\chi=0: 0.08%，: 0.08%，\chi=0.9: 1.5%），但 GPU 镜像存在未解决的 CPU↔WGSL 数值发散——渲染的盘没有引力透镜，尽管 CPU 追踪显示正确的 173° 偏折。

这是一个悬而未决的硬骨头：在 GPU 浮点精度下，精确 Kerr 哈密顿量的数值稳定性如何保证？两条等价实现为何发散？这或许是这个项目最深邃的开放问题。

结语：物理与工程的合奏

做这个渲染器的过程中，越来越确信一件事：物理保真与工程约束不是对立的取舍，而是可以互相成就的。

• 物理保真（自适应积分、Kerr 退化、ISCO、多普勒）通过 CPU 镜像变得可测试；
• 工程约束（GPU 不能挂起、WebGPU 规范、Bevy 0.19 的 bug）通过精心设计的语义（预算、底板、nudge）被驯服；
• 视觉目标（Gargantua）通过 HDR 管线、ACES、体积盘、抗混叠从物理中自然涌现，而非贴图伪造。

每一个像素都是一次从相机出发、穿越弯曲时空的微型宇宙学实验。六万亿次这样的实验每秒在 GPU 上发生。这就是实时黑洞渲染的诗意所在——用工程的严谨，让广义相对论的方程在眼前活过来。

附录：关键公式速查

本文基于 singularity-rs v0.2.0 源码（src/physics.rs、assets/shaders/black_hole.wgsl、src/render/）撰写。所有代码引用可在仓库中定位。