feat: Kerr deriv() with frame-dragging + spin-dependent horizon

assets/shaders/black_hole.wgsl

@@ -110,12 +110,21 @@ struct Deriv {
 fn deriv(pos: vec3<f32>, dir: vec3<f32>) -> Deriv {
     let r = length(pos);
     let rs = uniforms.rs;
+    // Kerr spin. χ ∈ [0,1]; a = χ·M, M = Rs/2 = 0.5 (Rs=1).
+    let chi = uniforms.spin;
+    let m = 0.5;
+    let a = chi * m;
+    // Schwarzschild radial bending (identical to Phase 1 at χ=0).
     let h = cross(pos, dir);
     let h2 = dot(h, h);
     let r5 = max(r * r * r * r * r, 1e-6);
-    let dpos = dir;
+    let radial = -1.5 * rs * h2 / r5 * pos;
-    let accel = -1.5 * rs * h2 / r5 * pos;
+    // Frame-dragging (Lense-Thirring leading term). Spin axis = +Y.
-    return Deriv(dpos, accel);
+    let spin_axis = vec3<f32>(0.0, 1.0, 0.0);
+    let r3 = max(r * r * r, 1e-6);
+    let drag = 2.0 * m * a / r3 * cross(spin_axis, dir);
+    let accel = radial + drag;
+    return Deriv(dir, accel);
 }
 
 // --- disk ---
@@ -267,7 +276,12 @@ fn fragment(in: VertexOutput) -> @location(0) vec4<f32> {
     var prev = pos;
     for (var i: u32 = 0u; i < steps; i = i + 1u) {
         let r = length(pos);
-        if (r < uniforms.rs) {
+        // Kerr horizon r+ = M + sqrt(M² - a²), M=0.5, a=χ·M. Equals Rs at χ=0.
+        let chi = uniforms.spin;
+        let m = 0.5;
+        let a = chi * m;
+        let r_plus = m + sqrt(max(m * m - a * a, 0.0));
+        if (r < r_plus) {
             // Captured: whatever we've composited so far is the result.
             break;
         }

src/physics.rs

@@ -72,6 +72,22 @@ pub fn kerr_horizon(chi: f32) -> f32 {
     m + (m * m - a * a).max(0.0).sqrt()
 }
 
+/// Kerr bending acceleration (CPU mirror of the shader `deriv` accel).
+/// `chi = a/M ∈ [0,1]`. At chi=0 this equals `bending_accel`.
+pub fn kerr_bending_accel(pos: Vec3, dir: Vec3, chi: f32) -> Vec3 {
+    let r = pos.length();
+    let m = 0.5;
+    let a = chi * m;
+    let h = pos.cross(dir);
+    let h2 = h.dot(h);
+    let r5 = (r * r * r * r * r).max(1e-6);
+    let radial = -1.5 * RS * h2 / r5 * pos;
+    let spin_axis = Vec3::Y;
+    let r3 = (r * r * r).max(1e-6);
+    let drag = 2.0 * m * a / r3 * spin_axis.cross(dir);
+    radial + drag
+}
+
 // Silence unused-import warning for Vec4 if not used; kept for future expansion.
 #[allow(dead_code)]
 fn _phantom(_v: Vec4) {}

tests/physics_test.rs

@@ -52,3 +52,25 @@ fn kerr_horizon_is_monotonically_decreasing() {
     assert!(a > b, "0.3 > 0.6: {} vs {}", a, b);
     assert!(b > c, "0.6 > 0.9: {} vs {}", b, c);
 }
+
+#[test]
+fn kerr_bending_accel_degenerates_to_schwarzschild_at_zero_spin() {
+    // At χ=0 the Kerr bending accel must equal the Schwarzschild one.
+    let pos = bevy::math::Vec3::new(3.0, 1.0, 4.0);
+    let dir = bevy::math::Vec3::new(0.2, -0.1, -0.97).normalize();
+    let schw = physics::bending_accel(pos, dir);
+    let kerr = physics::kerr_bending_accel(pos, dir, 0.0);
+    let diff = (schw - kerr).length();
+    assert!(diff < 1e-6, "spin=0 Kerr should match Schwarzschild; diff = {}", diff);
+}
+
+#[test]
+fn kerr_bending_accel_nonzero_off_axis_at_nonzero_spin() {
+    // At χ>0 the drag term must produce a different accel (frame-dragging exists).
+    let pos = bevy::math::Vec3::new(3.0, 1.0, 4.0);
+    let dir = bevy::math::Vec3::new(0.2, -0.1, -0.97).normalize();
+    let schw = physics::bending_accel(pos, dir);
+    let kerr = physics::kerr_bending_accel(pos, dir, 0.8);
+    let diff = (schw - kerr).length();
+    assert!(diff > 1e-4, "spin=0.8 Kerr should differ from Schwarzschild; diff = {}", diff);
+}