diff --git a/assets/shaders/black_hole.wgsl b/assets/shaders/black_hole.wgsl index 0be4aa3..16ab6ce 100644 --- a/assets/shaders/black_hole.wgsl +++ b/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, @@ -37,13 +38,59 @@ fn fragment(in: VertexOutput) -> @location(0) vec4 { 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(0.0, 0.0, 0.0, 1.0); + let steps = uniforms.steps; + + // Front-to-back compositing. + var accum_color = vec3(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(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(star, 1.0); + + return vec4(accum_color, 1.0); } diff --git a/assets/shaders/disk.wgsl b/assets/shaders/disk.wgsl new file mode 100644 index 0000000..330db3c --- /dev/null +++ b/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, a: f32) -> vec3 { + let c = cos(a); + let s = sin(a); + return vec3(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, cur: vec3) -> 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(t)); + let r = length(vec2(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, dir: vec3) -> vec3 { + let r = length(vec2(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(1.0, 0.95, 0.85), vec3(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(-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; +}