diff --git a/geodesic.comp b/geodesic.comp index cc4530f..9e7d1ad 100644 --- a/geodesic.comp +++ b/geodesic.comp @@ -30,15 +30,17 @@ layout(std140, binding = 3) uniform Objects { const float SagA_rs = 1.269e10; const float D_LAMBDA = 1e7; const double ESCAPE_R = 1e30; -vec4 objectColor = vec4(0.0); -int WIDTH = 400; -int HEIGHT = 300; +// Globals to store hit info +vec4 objectColor = vec4(0.0); +vec3 hitCenter = vec3(0.0); +float hitRadius = 0.0; struct Ray { float x, y, z, r, theta, phi; float dr, dtheta, dphi; - float E, L; }; + float E, L; +}; Ray initRay(vec3 pos, vec3 dir) { Ray ray; ray.x = pos.x; ray.y = pos.y; ray.z = pos.z; @@ -56,17 +58,22 @@ Ray initRay(vec3 pos, vec3 dir) { float dt_dL = sqrt((ray.dr*ray.dr)/f + ray.r*ray.r*(ray.dtheta*ray.dtheta + sin(ray.theta)*sin(ray.theta)*ray.dphi*ray.dphi)); ray.E = f * dt_dL; - return ray;} + return ray; +} bool intercept(Ray ray, float rs) { - return ray.r <= rs;} + return ray.r <= rs; +} +// Returns true on hit, captures center, radius, and base color bool interceptObject(Ray ray) { + vec3 P = vec3(ray.x, ray.y, ray.z); for (int i = 0; i < numObjects; ++i) { - vec3 objPos = objPosRadius[i].xyz; - float objRadius = objPosRadius[i].w; - float dist = distance(vec3(ray.x, ray.y, ray.z), objPos); - if (dist <= objRadius) { + vec3 center = objPosRadius[i].xyz; + float radius = objPosRadius[i].w; + if (distance(P, center) <= radius) { objectColor = objColor[i]; + hitCenter = center; + hitRadius = radius; return true; } } @@ -84,35 +91,12 @@ void geodesicRHS(Ray ray, out vec3 d1, out vec3 d2) { + (SagA_rs / (2.0 * r*r * f)) * dr * dr + r * (dtheta*dtheta + sin(theta)*sin(theta)*dphi*dphi); d2.y = -2.0*dr*dtheta/r + sin(theta)*cos(theta)*dphi*dphi; - d2.z = -2.0*dr*dphi/r - 2.0*cos(theta)/(sin(theta)) * dtheta * dphi;} + d2.z = -2.0*dr*dphi/r - 2.0*cos(theta)/(sin(theta)) * dtheta * dphi; +} void rk4Step(inout Ray ray, float dL) { - vec3 k1a, k1b, k2a, k2b, k3a, k3b, k4a, k4b; + vec3 k1a, k1b; geodesicRHS(ray, k1a, k1b); - //Ray r2; - //r2.r = ray.r + dL*0.5*k1a.x; r2.theta = ray.theta + dL*0.5*k1a.y; r2.phi = ray.phi + dL*0.5*k1a.z; - //r2.dr = ray.dr + dL*0.5*k1b.x; r2.dtheta = ray.dtheta + dL*0.5*k1b.y; r2.dphi = ray.dphi + dL*0.5*k1b.z; - //r2.E = ray.E; - //geodesicRHS(r2, k2a, k2b); -// - //Ray r3; - //r3.r = ray.r + dL*0.5*k2a.x; r3.theta = ray.theta + dL*0.5*k2a.y; r3.phi = ray.phi + dL*0.5*k2a.z; - //r3.dr = ray.dr + dL*0.5*k2b.x; r3.dtheta = ray.dtheta + dL*0.5*k2b.y; r3.dphi = ray.dphi + dL*0.5*k2b.z; - //r3.E = ray.E; - //geodesicRHS(r3, k3a, k3b); -// - //Ray r4; - //r4.r = ray.r + dL*k3a.x; r4.theta = ray.theta + dL*k3a.y; r4.phi = ray.phi + dL*k3a.z; - //r4.dr = ray.dr + dL*k3b.x; r4.dtheta = ray.dtheta + dL*k3b.y; r4.dphi = ray.dphi + dL*k3b.z; - //r4.E = ray.E; - //geodesicRHS(r4, k4a, k4b); -// - //ray.r += dL/6.0 * (k1a.x + 2*k2a.x + 2*k3a.x + k4a.x); - //ray.theta += dL/6.0 * (k1a.y + 2*k2a.y + 2*k3a.y + k4a.y); - //ray.phi += dL/6.0 * (k1a.z + 2*k2a.z + 2*k3a.z + k4a.z); - //ray.dr += dL/6.0 * (k1b.x + 2*k2b.x + 2*k3b.x + k4b.x); - //ray.dtheta += dL/6.0 * (k1b.y + 2*k2b.y + 2*k3b.y + k4b.y); - //ray.dphi += dL/6.0 * (k1b.z + 2*k2b.z + 2*k3b.z + k4b.z); ray.r += dL * k1a.x; ray.theta += dL * k1a.y; ray.phi += dL * k1a.z; @@ -122,25 +106,22 @@ void rk4Step(inout Ray ray, float dL) { ray.x = ray.r * sin(ray.theta) * cos(ray.phi); ray.y = ray.r * sin(ray.theta) * sin(ray.phi); - ray.z = ray.r * cos(ray.theta);} + ray.z = ray.r * cos(ray.theta); +} bool crossesEquatorialPlane(vec3 oldPos, vec3 newPos) { bool crossed = (oldPos.y * newPos.y < 0.0); float r = length(vec2(newPos.x, newPos.z)); - return crossed && (r >= disk_r1 && r <= disk_r2);} - + return crossed && (r >= disk_r1 && r <= disk_r2); +} void main() { - //if (cam.moving) { - // WIDTH = 200; - // HEIGHT = 150; - //} else { - // WIDTH = 800; - // HEIGHT = 600; - //} + int WIDTH = cam.moving ? 200 : 200; + int HEIGHT = cam.moving ? 150 : 150; + ivec2 pix = ivec2(gl_GlobalInvocationID.xy); if (pix.x >= WIDTH || pix.y >= HEIGHT) return; - // -- Init Ray -- // + // Init Ray float u = (2.0 * (pix.x + 0.5) / WIDTH - 1.0) * cam.aspect * cam.tanHalfFov; float v = (1.0 - 2.0 * (pix.y + 0.5) / HEIGHT) * cam.tanHalfFov; vec3 dir = normalize(u * cam.camRight - v * cam.camUp + cam.camForward); @@ -149,52 +130,45 @@ void main() { vec4 color = vec4(0.0); vec3 prevPos = vec3(ray.x, ray.y, ray.z); float lambda = 0.0; - // -- Intercepts -- // + bool hitBlackHole = false; bool hitDisk = false; bool hitObject = false; - int steps; - if (cam.moving) { - steps=20000; - } else { - steps=80000; - }; + int steps = cam.moving ? 60000 : 60000; - // Step Loop for (int i = 0; i < steps; ++i) { - if (intercept(ray, SagA_rs)) { - hitBlackHole = true; - break; - } + if (intercept(ray, SagA_rs)) { hitBlackHole = true; break; } rk4Step(ray, D_LAMBDA); - lambda += D_LAMBDA; - vec3 newPos = vec3(ray.x, ray.y, ray.z); - if (crossesEquatorialPlane(prevPos, newPos)) { - float r = length(vec2(newPos.x, newPos.z)); - hitDisk = true; - break; - } - - if (interceptObject(ray)) { - hitObject = true; - break; - } + vec3 newPos = vec3(ray.x, ray.y, ray.z); + if (crossesEquatorialPlane(prevPos, newPos)) { hitDisk = true; break; } + if (interceptObject(ray)) { hitObject = true; break; } prevPos = newPos; if (ray.r > ESCAPE_R) break; } if (hitDisk) { - double r = sqrt(ray.x*ray.x + ray.y*ray.y + ray.z*ray.z); - r = r / (SagA_rs * 4.2); - vec3 diskColor = vec3(1.0f, r, 0.2); - color = vec4(diskColor * 1.0, 1.0); + double r = length(vec3(ray.x, ray.y, ray.z)) / disk_r2; + vec3 diskColor = vec3(1.0, r, 0.2); + //r = 1.0 - abs(r - 0.5) * 2.0; + color = vec4(diskColor, r); + } else if (hitBlackHole) { color = vec4(0.0, 0.0, 0.0, 1.0); - } else if (hitObject){ - color = objectColor; + + } else if (hitObject) { + // Compute shading + vec3 P = vec3(ray.x, ray.y, ray.z); + vec3 N = normalize(P - hitCenter); + vec3 V = normalize(cam.camPos - P); + float ambient = 0.1; + float diff = max(dot(N, V), 0.0); + float intensity = ambient + (1.0 - ambient) * diff; + vec3 shaded = objectColor.rgb * intensity; + color = vec4(shaded, objectColor.a); + } else { color = vec4(0.0); }