feat: accretion disk with Doppler beaming + lensed Einstein halo

assets/shaders/black_hole.wgsl

@@ -2,6 +2,7 @@
 #import "shaders/ray_gen.wgsl"
 #import "shaders/geodesic_schwarzschild.wgsl"
 #import "shaders/stars.wgsl"
+#import "shaders/disk.wgsl"
 
 struct BlackHoleUniforms {
     eye: vec4<f32>,
@@ -37,13 +38,59 @@ fn fragment(in: VertexOutput) -> @location(0) vec4<f32> {
     var uv = (in.uv * 2.0 - 1.0);
     uv.x *= aspect;
     let dir = ray_direction(uv);
+
+    // Work in disk-local space: rotate eye + dir by -disk_tilt around X so the
+    // disk lies on y=0. (disk_hit/disk_color assume disk-local coords.)
+    var pos = rot_x(uniforms.eye.xyz, -uniforms.disk_tilt);
+    var d   = normalize(rot_x(dir, -uniforms.disk_tilt));
+
     let dt = max(length(uniforms.eye.xyz), 20.0) / f32(uniforms.steps);
-    let res = classify_ray(uniforms.eye.xyz, dir, uniforms.steps, dt);
-    if (res.status == 1u) {
-        // Captured = shadow.
-        return vec4<f32>(0.0, 0.0, 0.0, 1.0);
+    let steps = uniforms.steps;
+
+    // Front-to-back compositing.
+    var accum_color = vec3<f32>(0.0);
+    var accum_alpha = 0.0;
+
+    var prev = pos;
+    for (var i: u32 = 0u; i < steps; i = i + 1u) {
+        let r = length(pos);
+        if (r < uniforms.rs) {
+            // Captured: whatever we've composited so far is the result.
+            break;
+        }
+        if (r > 1000.0) {
+            // Escaped: add background stars along the (disk-local) final dir.
+            // Rotate back to world for the star sample.
+            let world_dir = normalize(rot_x(d, uniforms.disk_tilt));
+            let star = star_color(world_dir, uniforms.star_intensity);
+            accum_color += (1.0 - accum_alpha) * star;
+            accum_alpha = 1.0;
+            break;
+        }
+
+        // RK4 step (single step), then test disk crossing on the segment.
+        let (k1p, k1d) = deriv(pos, d);
+        let (k2p, k2d) = deriv(pos + k1p * dt * 0.5, normalize(d + k1d * dt * 0.5));
+        let (k3p, k3d) = deriv(pos + k2p * dt * 0.5, normalize(d + k2d * dt * 0.5));
+        let (k4p, k4d) = deriv(pos + k3p * dt, normalize(d + k3d * dt));
+        let new_pos = pos + (k1p + 2.0*k2p + 2.0*k3p + k4p) * dt / 6.0;
+        let new_dir = normalize(d + (k1d + 2.0*k2d + 2.0*k3d + k4d) * dt / 6.0);
+
+        if (disk_hit(prev, new_pos)) {
+            // Approximate the crossing point by interpolating to y=0.
+            let ty = prev.y / (prev.y - new_pos.y);
+            let hit = mix(prev, new_pos, vec3<f32>(ty));
+            let dc = disk_color(hit, new_dir);
+            let a = 0.85; // disk is nearly opaque
+            accum_color += (1.0 - accum_alpha) * dc * a;
+            accum_alpha += (1.0 - accum_alpha) * a;
+            if (accum_alpha > 0.99) { break; }
+        }
+
+        prev = new_pos;
+        pos = new_pos;
+        d = new_dir;
     }
-    // Escaped: sample procedural stars along the bent final direction.
-    let star = star_color(normalize(res.final_dir), uniforms.star_intensity);
-    return vec4<f32>(star, 1.0);
+
+    return vec4<f32>(accum_color, 1.0);
 }

assets/shaders/disk.wgsl

@@ -0,0 +1,61 @@
+// Disk plane is the xz-plane in world space, tilted by `disk_tilt` around the
+// x-axis. We work in "disk-local" coordinates by rotating the ray.
+
+// Rotate a vector around the X axis by angle a.
+fn rot_x(v: vec3<f32>, a: f32) -> vec3<f32> {
+    let c = cos(a);
+    let s = sin(a);
+    return vec3<f32>(v.x, c * v.y - s * v.z, s * v.y + c * v.z);
+}
+
+// Returns true if the segment pos->pos+dir*dt crosses the disk plane (y=0)
+// within radius [disk_inner, disk_outer]. (prev, cur are the segment endpoints.)
+fn disk_hit(prev: vec3<f32>, cur: vec3<f32>) -> bool {
+    let y0 = prev.y;
+    let y1 = cur.y;
+    if (y0 * y1 > 0.0) {
+        return false; // same side, no crossing
+    }
+    // Linear interpolate to the crossing point.
+    let t = y0 / (y0 - y1);
+    let cross = mix(prev, cur, vec3<f32>(t));
+    let r = length(vec2<f32>(cross.x, cross.z));
+    return r >= uniforms.disk_inner && r <= uniforms.disk_outer;
+}
+
+// Shade a disk hit: procedural texture + Doppler beaming + temperature color.
+fn disk_color(pos: vec3<f32>, dir: vec3<f32>) -> vec3<f32> {
+    let r = length(vec2<f32>(pos.x, pos.z));
+    let phi = atan2(pos.z, pos.x);
+
+    // Procedural noise: layered angular + radial, animated by rotation.
+    let rot = uniforms.time * uniforms.disk_rotation_speed / pow(r, 1.5);
+    let n = sin(phi * 8.0 + rot) * 0.5 + 0.5;
+    let n2 = sin(phi * 23.0 - rot * 1.7 + r * 2.0) * 0.5 + 0.5;
+    let noise = mix(n, n2, 0.4);
+
+    // Temperature gradient: hotter (white-blue) near inner edge, cooler (orange-red) outer.
+    let t = (r - uniforms.disk_inner) / (uniforms.disk_outer - uniforms.disk_inner);
+    let tcol = mix(vec3<f32>(1.0, 0.95, 0.85), vec3<f32>(1.0, 0.45, 0.12), clamp(t, 0.0, 1.0));
+
+    // Falloff: brighter at inner edge.
+    let falloff = 1.0 / pow(r / uniforms.disk_inner, 2.0);
+
+    var col = tcol * (0.6 + 0.4 * noise) * falloff;
+
+    // Doppler beaming. Disk orbits Keplerian-ish: v ~ sqrt(Rs/(2r)).
+    let v_orbital = sqrt(uniforms.rs / (2.0 * r));
+    // Orbital velocity direction (tangent) in the disk plane.
+    let tangent = normalize(vec3<f32>(-sin(phi), 0.0, cos(phi)));
+    // Scalar approximation: projection of orbital velocity onto ray direction.
+    let vdotn = dot(tangent * v_orbital, -dir); // toward viewer if positive
+    let gamma = 1.0 / sqrt(max(1.0 - v_orbital * v_orbital, 1e-4));
+    var doppler = 1.0;
+    if (uniforms.doppler_enabled != 0u) {
+        let delta = 1.0 / (gamma * (1.0 - vdotn));
+        doppler = pow(delta, 3.0) * uniforms.doppler_strength;
+    }
+    col *= doppler;
+
+    return col * uniforms.disk_brightness;
+}