diff --git a/assets/shaders/black_hole.wgsl b/assets/shaders/black_hole.wgsl index 455e4c6..d74f2b3 100644 --- a/assets/shaders/black_hole.wgsl +++ b/assets/shaders/black_hole.wgsl @@ -110,12 +110,21 @@ struct Deriv { fn deriv(pos: vec3, dir: vec3) -> 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 accel = -1.5 * rs * h2 / r5 * pos; - return Deriv(dpos, accel); + let radial = -1.5 * rs * h2 / r5 * pos; + // Frame-dragging (Lense-Thirring leading term). Spin axis = +Y. + let spin_axis = vec3(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 { 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; } diff --git a/src/physics.rs b/src/physics.rs index 6751e34..71d30f2 100644 --- a/src/physics.rs +++ b/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) {} diff --git a/tests/physics_test.rs b/tests/physics_test.rs index 45e0fc1..5018793 100644 --- a/tests/physics_test.rs +++ b/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); +}