diff --git a/2D_lensing.cpp b/2D_lensing.cpp deleted file mode 100644 index f2f4e79..0000000 --- a/2D_lensing.cpp +++ /dev/null @@ -1,231 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#define _USE_MATH_DEFINES -#include -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -using namespace glm; -using namespace std; - -double c = 299792458.0; -double G = 6.67430e-11; - -struct Ray; -void rk4Step(Ray& ray, double dλ, double rs); - -// --- Structs --- // -struct Engine { - GLFWwindow* window; - int WIDTH = 800; - int HEIGHT = 600; - float width = 100000000000.0f; // Width of the viewport in meters - float height = 75000000000.0f; // Height of the viewport in meters - - // Navigation state - float offsetX = 0.0f, offsetY = 0.0f; - float zoom = 1.0f; - bool middleMousePressed = false; - double lastMouseX = 0.0, lastMouseY = 0; - - Engine() { - if (!glfwInit()) { - cerr << "Failed to initialize GLFW" << endl; - exit(EXIT_FAILURE); - } - window = glfwCreateWindow(WIDTH, HEIGHT, "Black Hole Simulation", NULL, NULL); - if (!window) { - cerr << "Failed to create GLFW window" << endl; - glfwTerminate(); - exit(EXIT_FAILURE); - } - glfwMakeContextCurrent(window); - glewExperimental = GL_TRUE; - if (glewInit() != GLEW_OK) { - cerr << "Failed to initialize GLEW" << endl; - glfwDestroyWindow(window); - glfwTerminate(); - exit(EXIT_FAILURE); - } - glViewport(0, 0, WIDTH, HEIGHT);; - } - - void run() { - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - glMatrixMode(GL_PROJECTION); - glLoadIdentity(); - double left = -width + offsetX; - double right = width + offsetX; - double bottom = -height + offsetY; - double top = height + offsetY; - glOrtho(left, right, bottom, top, -1.0, 1.0); - glMatrixMode(GL_MODELVIEW); - glLoadIdentity(); - } -}; -Engine engine; -struct BlackHole { - vec3 position; - double mass; - double radius; - double r_s; - - BlackHole(vec3 pos, float m) : position(pos), mass(m) {r_s = 2.0 * G * mass / (c*c);} - void draw() { - glBegin(GL_TRIANGLE_FAN); - glColor3f(1.0f, 0.0f, 0.0f); // Red color for the black hole - glVertex2f(0.0f, 0.0f); // Center - for(int i = 0; i <= 100; i++) { - float angle = 2.0f * M_PI * i / 100; - float x = r_s * cos(angle); // Radius of 0.1 - float y = r_s * sin(angle); - glVertex2f(x, y); - } - glEnd(); - } -}; -BlackHole SagA(vec3(0.0f, 0.0f, 0.0f), 8.54e36); // Sagittarius A black hole -struct Ray{ - // -- cartesian coords -- // - double x; double y; - // -- polar coords -- // - double r; double phi; - double dr; double dphi; - vector trail; // trail of points - double E, L; // conserved quantities - - Ray(vec2 pos, vec2 dir) : x(pos.x), y(pos.y), r(sqrt(pos.x * pos.x + pos.y * pos.y)), phi(atan2(pos.y, pos.x)), dr(dir.x), dphi(dir.y) { - // step 1) get polar coords (r, phi) : - this->r = sqrt(x*x + y*y); - this->phi = atan2(y, x); - // step 2) seed velocities : - dr = dir.x * cos(phi) + dir.y * sin(phi); // m/s - dphi = ( -dir.x * sin(phi) + dir.y * cos(phi) ) / r; - // step 3) store conserved quantities - L = r*r * dphi; - double f = 1.0 - SagA.r_s/r; - double dt_dλ = sqrt( (dr*dr)/(f*f) + (r*r*dphi*dphi)/f ); - E = f * dt_dλ; - // step 4) start trail : - trail.push_back({x, y}); - } - void draw(const std::vector& rays) { - // draw current ray positions as points - glPointSize(2.0f); - glColor3f(1.0f, 0.0f, 0.0f); - glBegin(GL_POINTS); - for (const auto& ray : rays) { - glVertex2f(ray.x, ray.y); - } - glEnd(); - - // turn on blending for the trails - glEnable(GL_BLEND); - glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); - glLineWidth(1.0f); - - // draw each trail with fading alpha - for (const auto& ray : rays) { - size_t N = ray.trail.size(); - if (N < 2) continue; - - glBegin(GL_LINE_STRIP); - for (size_t i = 0; i < N; ++i) { - // older points (i=0) get alpha≈0, newer get alpha≈1 - float alpha = float(i) / float(N - 1); - glColor4f(1.0f, 1.0f, 1.0f, std::max(alpha, 0.05f)); - glVertex2f(ray.trail[i].x, ray.trail[i].y); - } - glEnd(); - } - - glDisable(GL_BLEND); - } - void step(double dλ, double rs) { - // 1) integrate (r,φ,dr,dφ) - if(r <= rs) return; // stop if inside the event horizon - rk4Step(*this, dλ, rs); - - // 2) convert back to cartesian x,y - x = r * cos(phi); - y = r * sin(phi); - - // 3) record the trail - trail.push_back({ float(x), float(y) }); - } -}; -vector rays; - -void geodesicRHS(const Ray& ray, double rhs[4], double rs) { - double r = ray.r; - double dr = ray.dr; - double dphi = ray.dphi; - double E = ray.E; - - double f = 1.0 - rs/r; - - // dr/dλ = dr - rhs[0] = dr; - // dφ/dλ = dphi - rhs[1] = dphi; - - // d²r/dλ² from Schwarzschild null geodesic: - double dt_dλ = E / f; - rhs[2] = - - (rs/(2*r*r)) * f * (dt_dλ*dt_dλ) - + (rs/(2*r*r*f)) * (dr*dr) - + (r - rs) * (dphi*dphi); - - // d²φ/dλ² = -2*(dr * dphi) / r - rhs[3] = -2.0 * dr * dphi / r; -} -void addState(const double a[4], const double b[4], double factor, double out[4]) { - for (int i = 0; i < 4; i++) - out[i] = a[i] + b[i] * factor; -} -void rk4Step(Ray& ray, double dλ, double rs) { - double y0[4] = { ray.r, ray.phi, ray.dr, ray.dphi }; - double k1[4], k2[4], k3[4], k4[4], temp[4]; - - geodesicRHS(ray, k1, rs); - addState(y0, k1, dλ/2.0, temp); - Ray r2 = ray; r2.r=temp[0]; r2.phi=temp[1]; r2.dr=temp[2]; r2.dphi=temp[3]; - geodesicRHS(r2, k2, rs); - - addState(y0, k2, dλ/2.0, temp); - Ray r3 = ray; r3.r=temp[0]; r3.phi=temp[1]; r3.dr=temp[2]; r3.dphi=temp[3]; - geodesicRHS(r3, k3, rs); - - addState(y0, k3, dλ, temp); - Ray r4 = ray; r4.r=temp[0]; r4.phi=temp[1]; r4.dr=temp[2]; r4.dphi=temp[3]; - geodesicRHS(r4, k4, rs); - - ray.r += (dλ/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); - ray.phi += (dλ/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]); - ray.dr += (dλ/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]); - ray.dphi += (dλ/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]); -} - - -int main () { - //rays.push_back(Ray(vec2(-1e11, 3.27606302719999999e10), vec2(c, 0.0f))); - while(!glfwWindowShouldClose(engine.window)) { - engine.run(); - SagA.draw(); - - for (auto& ray : rays) { - ray.step(1.0f, SagA.r_s); - ray.draw(rays); - } - - glfwSwapBuffers(engine.window); - glfwPollEvents(); - } - - return 0; -} diff --git a/CPU-geodesic.cpp b/CPU-geodesic.cpp deleted file mode 100644 index e1f7e30..0000000 --- a/CPU-geodesic.cpp +++ /dev/null @@ -1,496 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#define _USE_MATH_DEFINES -#include -#include -#include -#include -#include -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -using namespace glm; -using namespace std; -using Clock = std::chrono::high_resolution_clock; - -// VARS -double lastPrintTime = 0.0; -int framesCount = 0; -double c = 299792458.0; -double G = 6.67430e-11; -bool useGeodesics = false; - -struct Camera { - vec3 pos; - vec3 target; - float fovY; - float azimuth, elevation, radius; - float minRadius = 1e12f, maxRadius = 1e20f; - bool dragging = false; - bool panning = false; - double lastX = 0, lastY = 0; - - // Adjustable speeds - float orbitSpeed = 0.008f; - float panSpeed = 0.001f; - float zoomSpeed = 1.08f; // closer to 1 = slower zoom - - Camera() : azimuth(0), elevation(M_PI / 2.0f), radius(6.34194e10), fovY(60.0f) { - target = vec3(0, 0, 0); - updateVectors(); - } - - void updateVectors() { - pos.x = target.x + radius * sin(elevation) * cos(azimuth); - pos.y = target.y + radius * cos(elevation); - pos.z = target.z + radius * sin(elevation) * sin(azimuth); - } - void processMouse(GLFWwindow* window, double xpos, double ypos) { - float dx = float(xpos - lastX), dy = float(ypos - lastY); - if (dragging && !panning) { - // Orbit - azimuth -= dx * orbitSpeed; - elevation -= dy * orbitSpeed; - elevation = clamp(elevation, 0.01f, float(M_PI)-0.01f); - } else if (panning) { - // Pan (move target in camera plane) - vec3 forward = normalize(target - pos); - vec3 right = normalize(cross(forward, vec3(0,1,0))); - vec3 up = cross(right, forward); - target += -right * dx * panSpeed * radius + up * dy * panSpeed * radius; - } - updateVectors(); - lastX = xpos; lastY = ypos; - } - void processScroll(double yoffset) { - // Zoom (dolly in/out) - if (yoffset < 0) - radius *= pow(zoomSpeed, -yoffset); - else - radius /= pow(zoomSpeed, yoffset); - radius = clamp(radius, minRadius, maxRadius); - updateVectors(); - } - static void mouseButtonCallback(GLFWwindow* window, int button, int action, int mods) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(window); - if (button == GLFW_MOUSE_BUTTON_LEFT) { - if (action == GLFW_PRESS) { - cam->dragging = true; - cam->panning = (mods & GLFW_MOD_SHIFT); - double x, y; glfwGetCursorPos(window, &x, &y); - cam->lastX = x; cam->lastY = y; - } else if (action == GLFW_RELEASE) { - cam->dragging = false; - cam->panning = false; - } - } - } - static void cursorPosCallback(GLFWwindow* window, double xpos, double ypos) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(window); - cam->processMouse(window, xpos, ypos); - } - static void scrollCallback(GLFWwindow* window, double xoffset, double yoffset) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(window); - cam->processScroll(yoffset); - } -}; -Camera camera; - -struct Ray; -void rk4Step(Ray& ray, double dλ, double rs); - -struct Engine { - // -- Quad & Texture render -- // - GLFWwindow* window; - GLuint quadVAO; - GLuint texture; - GLuint shaderProgram; - int WIDTH = 800; - int HEIGHT = 600; - float width = 100000000000.0f; // Width of the viewport in meters - float height = 75000000000.0f; // Height of the viewport in meters - - Engine() { - if (!glfwInit()) { - cerr << "GLFW init failed\n"; - exit(EXIT_FAILURE); - } - glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3); - glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3); - glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE); - window = glfwCreateWindow(WIDTH, HEIGHT, "Black Hole", nullptr, nullptr); - if (!window) { - cerr << "Failed to create GLFW window\n"; - glfwTerminate(); - exit(EXIT_FAILURE); - } - glfwMakeContextCurrent(window); - glewExperimental = GL_TRUE; - GLenum glewErr = glewInit(); - if (glewErr != GLEW_OK) { - cerr << "Failed to initialize GLEW: " - << (const char*)glewGetErrorString(glewErr) - << "\n"; - glfwTerminate(); - exit(EXIT_FAILURE); - } - cout << "OpenGL " << glGetString(GL_VERSION) << "\n"; - this->shaderProgram = CreateShaderProgram(); - - auto result = QuadVAO(); - this->quadVAO = result[0]; - this->texture = result[1]; - } - GLuint CreateShaderProgram(){ - const char* vertexShaderSource = R"( - #version 330 core - layout (location = 0) in vec2 aPos; // Changed to vec2 - layout (location = 1) in vec2 aTexCoord; - out vec2 TexCoord; - void main() { - gl_Position = vec4(aPos, 0.0, 1.0); // Explicit z=0 - TexCoord = aTexCoord; - })"; - - const char* fragmentShaderSource = R"( - #version 330 core - in vec2 TexCoord; - out vec4 FragColor; - uniform sampler2D screenTexture; - void main() { - FragColor = texture(screenTexture, TexCoord); - })"; - - // vertex shader - GLuint vertexShader = glCreateShader(GL_VERTEX_SHADER); - glShaderSource(vertexShader, 1, &vertexShaderSource, nullptr); - glCompileShader(vertexShader); - - // fragment shader - GLuint fragmentShader = glCreateShader(GL_FRAGMENT_SHADER); - glShaderSource(fragmentShader, 1, &fragmentShaderSource, nullptr); - glCompileShader(fragmentShader); - - GLuint shaderProgram = glCreateProgram(); - glAttachShader(shaderProgram, vertexShader); - glAttachShader(shaderProgram, fragmentShader); - glLinkProgram(shaderProgram); - - glDeleteShader(vertexShader); - glDeleteShader(fragmentShader); - - return shaderProgram; - }; - vector QuadVAO(){ - float quadVertices[] = { - // positions // texCoords - -1.0f, 1.0f, 0.0f, 1.0f, // top left - -1.0f, -1.0f, 0.0f, 0.0f, // bottom left - 1.0f, -1.0f, 1.0f, 0.0f, // bottom right - - -1.0f, 1.0f, 0.0f, 1.0f, // top left - 1.0f, -1.0f, 1.0f, 0.0f, // bottom right - 1.0f, 1.0f, 1.0f, 1.0f // top right - - }; - - GLuint VAO, VBO; - glGenVertexArrays(1, &VAO); - glGenBuffers(1, &VBO); - - glBindVertexArray(VAO); - glBindBuffer(GL_ARRAY_BUFFER, VBO); - glBufferData(GL_ARRAY_BUFFER, sizeof(quadVertices), quadVertices, GL_STATIC_DRAW); - - glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)0); - glEnableVertexAttribArray(0); - glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)(2 * sizeof(float))); - glEnableVertexAttribArray(1); - - GLuint texture; - glGenTextures(1, &texture); - glBindTexture(GL_TEXTURE_2D, texture); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - vector VAOtexture = {VAO, texture}; - return VAOtexture; - } - void renderScene(const vector& pixels, int texWidth, int texHeight) { - // update texture w/ ray-tracing results - glBindTexture(GL_TEXTURE_2D, texture); - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, texWidth, texHeight, 0, GL_RGB, GL_UNSIGNED_BYTE, pixels.data()); - - // clear screen and draw textured quad - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - glUseProgram(shaderProgram); - - GLint textureLocation = glGetUniformLocation(shaderProgram, "screenTexture"); - glUniform1i(textureLocation, 0); - - glBindVertexArray(quadVAO); - glDrawArrays(GL_TRIANGLES, 0, 6); - - glfwSwapBuffers(window); - glfwPollEvents(); - }; - static void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods) { - if (action == GLFW_PRESS) { - if (key == GLFW_KEY_G) { - useGeodesics = !useGeodesics; - cout << "Geodesics: " << (useGeodesics ? "ON\n" : "OFF\n"); - } - } - } -}; -Engine engine; -struct BlackHole { - vec3 position; - double mass; - double radius; - double r_s; - - BlackHole(vec3 pos, float m) : position(pos), mass(m) {r_s = 2.0 * G * mass / (c*c);} - bool Intercept(float px, float py, float pz) const { - float dx = px - position.x; - float dy = py - position.y; - float dz = pz - position.z; - float dist2 = dx * dx + dy * dy + dz * dz; - return dist2 < r_s * r_s; - } -}; -BlackHole SagA(vec3(0.0f, 0.0f, 0.0f), 8.54e36); // Sagittarius A black hole -struct Ray{ - // -- cartesian coords -- // - double x; double y; double z; - // -- polar coords -- // - double r; double phi; double theta; - double dr; double dphi; double dtheta; - double E, L; // conserved quantities - - Ray(vec3 pos, vec3 dir) : x(pos.x), y(pos.y), z(pos.z) { - // Step 1: get spherical coords (r, theta, phi) - r = sqrt(x*x + y*y + z*z); - theta = acos(z / r); - phi = atan2(y, x); - - // Step 2: seed velocities (dr, dtheta, dphi) - // Convert direction to spherical basis - double dx = dir.x, dy = dir.y, dz = dir.z; - dr = sin(theta)*cos(phi)*dx + sin(theta)*sin(phi)*dy + cos(theta)*dz; - dtheta = cos(theta)*cos(phi)*dx + cos(theta)*sin(phi)*dy - sin(theta)*dz; - dtheta /= r; - dphi = -sin(phi)*dx + cos(phi)*dy; - dphi /= (r * sin(theta)); - - // Step 3: store conserved quantities - L = r * r * sin(theta) * dphi; - double f = 1.0 - SagA.r_s / r; - double dt_dλ = sqrt((dr*dr)/f + r*r*dtheta*dtheta + r*r*sin(theta)*sin(theta)*dphi*dphi); - E = f * dt_dλ; - } - void step(double dλ, double rs) { - if (r <= rs) return; - rk4Step(*this, dλ, rs); - // convert back to cartesian - this->x = r * sin(theta) * cos(phi); - this->y = r * sin(theta) * sin(phi); - this->z = r * cos(theta); - } -}; - -void raytrace(vector& pixels, int W, int H) { - pixels.resize(W * H * 3); - - // build camera basis - vec3 forward = normalize(camera.target - camera.pos); - vec3 right = normalize(cross(forward, vec3(0,1,0))); - vec3 up = cross(right, forward); - float aspect = float(W) / float(H); - float tanHalfFov = tan(radians(camera.fovY) * 0.5f); - - #pragma omp parallel for schedule(dynamic, 4) - for(int y = 0; y < H; ++y) { - for(int x = 0; x < W; ++x) { - // NDC → screen space in [−1,1] - float u = (2.0f * (x + 0.5f) / float(W) - 1.0f) * aspect * tanHalfFov; - float v = (1.0f - 2.0f * (y + 0.5f) / float(H)) * tanHalfFov; - vec3 dir = normalize(u*right + v*up + forward); - - // construct your Ray - Ray ray(camera.pos, dir); - - const int MAX_STEPS = 10000; - const double D_LAMBDA = 1e7; - const double ESCAPE_R = 1e14; - - // 2) march the ray forward in λ - vec3 color(0.0f); - if (!useGeodesics) { - double b = 2.0 * dot(camera.pos, dir); - double c0 = dot(camera.pos, camera.pos) - SagA.r_s*SagA.r_s; - double disc = b*b - 4.0*c0; - if (disc > 0.0) { - double t1 = (-b - sqrt(disc)) * 0.5; - double t2 = (-b + sqrt(disc)) * 0.5; - if (t1 > 0.0 || t2 > 0.0) - color = vec3(1.0f, 0.0f, 0.0f); - } - } - else { - // full null‐geodesic march - Ray ray(camera.pos, dir); - for(int i = 0; i < MAX_STEPS; ++i) { - if (SagA.Intercept(ray.x, ray.y, ray.z)) { - color = vec3(1.0f, 0.0f, 0.0f); - break; - } - ray.step(D_LAMBDA, SagA.r_s); - if (ray.r > ESCAPE_R) { - // escaped to infinity → remains black - break; - } - } - } - - int idx = (y * W + x) * 3; - pixels[idx+0] = (unsigned char)(color.r * 255); - pixels[idx+1] = (unsigned char)(color.g * 255); - pixels[idx+2] = (unsigned char)(color.b * 255); - } - } -} - -void geodesicRHS(const Ray& ray, double rhs[6], double rs) { - double r = ray.r; - double theta = ray.theta; - double dr = ray.dr; - double dtheta = ray.dtheta; - double dphi = ray.dphi; - double E = ray.E; - - double f = 1.0 - rs / r; - double dt_dlambda = E / f; - - // First derivatives - rhs[0] = dr; - rhs[1] = dtheta; - rhs[2] = dphi; - - // Second derivatives (from 3D Schwarzschild null geodesics): - rhs[3] = - - (rs / (2 * r * r)) * f * dt_dlambda * dt_dlambda - + (rs / (2 * r * r * f)) * dr * dr - + r * (dtheta * dtheta + sin(theta) * sin(theta) * dphi * dphi); - - rhs[4] = - - (2.0 / r) * dr * dtheta - + sin(theta) * cos(theta) * dphi * dphi; - - rhs[5] = - - (2.0 / r) * dr * dphi - - 2.0 * cos(theta) / sin(theta) * dtheta * dphi; -} -void addState(const double a[6], const double b[6], double factor, double out[6]) { - for (int i = 0; i < 6; i++) - out[i] = a[i] + b[i] * factor; -} -void rk4Step(Ray& ray, double dλ, double rs) { - double y0[6] = { ray.r, ray.theta, ray.phi, ray.dr, ray.dtheta, ray.dphi }; - double k1[6], k2[6], k3[6], k4[6], temp[6]; - - geodesicRHS(ray, k1, rs); - addState(y0, k1, dλ/2.0, temp); - Ray r2 = ray; - r2.r = temp[0]; r2.theta = temp[1]; r2.phi = temp[2]; - r2.dr = temp[3]; r2.dtheta = temp[4]; r2.dphi = temp[5]; - geodesicRHS(r2, k2, rs); - - addState(y0, k2, dλ/2.0, temp); - Ray r3 = ray; - r3.r = temp[0]; r3.theta = temp[1]; r3.phi = temp[2]; - r3.dr = temp[3]; r3.dtheta = temp[4]; r3.dphi = temp[5]; - geodesicRHS(r3, k3, rs); - - addState(y0, k3, dλ, temp); - Ray r4 = ray; - r4.r = temp[0]; r4.theta = temp[1]; r4.phi = temp[2]; - r4.dr = temp[3]; r4.dtheta = temp[4]; r4.dphi = temp[5]; - geodesicRHS(r4, k4, rs); - - ray.r += (dλ/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); - ray.theta += (dλ/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]); - ray.phi += (dλ/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]); - ray.dr += (dλ/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]); - ray.dtheta += (dλ/6.0)*(k1[4] + 2*k2[4] + 2*k3[4] + k4[4]); - ray.dphi += (dλ/6.0)*(k1[5] + 2*k2[5] + 2*k3[5] + k4[5]); -} - -void setupCameraCallbacks(GLFWwindow* window) { - glfwSetWindowUserPointer(window, &camera); - glfwSetMouseButtonCallback(window, Camera::mouseButtonCallback); - glfwSetCursorPosCallback(window, Camera::cursorPosCallback); - glfwSetScrollCallback(window, Camera::scrollCallback); - glfwSetKeyCallback(window, Engine::keyCallback); -} - -// -- MAIN -- // -int main() { - setupCameraCallbacks(engine.window); - vector pixels(engine.WIDTH * engine.HEIGHT * 3); - - auto t0 = Clock::now(); - lastPrintTime = std::chrono::duration(t0.time_since_epoch()).count(); - - while (!glfwWindowShouldClose(engine.window)) { - raytrace(pixels, engine.WIDTH, engine.HEIGHT); - engine.renderScene(pixels, engine.WIDTH, engine.HEIGHT); - - // 2) FPS counting - framesCount++; - auto t1 = Clock::now(); - double now = std::chrono::duration(t1.time_since_epoch()).count(); - if (now - lastPrintTime >= 1.0) { - cout << "FPS: " << framesCount / (now - lastPrintTime) << "\n"; - framesCount = 0; - lastPrintTime = now; - } - - } - - glfwDestroyWindow(engine.window); - glfwTerminate(); - return 0; -} - - - - - - - - - - - - - - - - - - - // 2) FPS counting - // framesCount++; - // auto t1 = Clock::now(); - // double now = std::chrono::duration(t1.time_since_epoch()).count(); - // if (now - lastPrintTime >= 1.0) { - // cout << "FPS: " << framesCount / (now - lastPrintTime) << "\n"; - // framesCount = 0; - // lastPrintTime = now; - // } - //raytrace(pixels, engine.WIDTH, engine.HEIGHT); diff --git a/Gravity_Sim.zip b/Gravity_Sim.zip deleted file mode 100644 index 098c42d..0000000 Binary files a/Gravity_Sim.zip and /dev/null differ diff --git a/README.md b/README.md index aae4faf..760b012 100644 --- a/README.md +++ b/README.md @@ -1,29 +1,2 @@ # black_hole -Black hole simulation project -Here is the black hole raw code, everything will be inside a src bin incase you want to copy the files -I'm writing this as I'm beginning this project (hopefully I complete it ;D) here is what I plan to do: - -1. Ray-tracing : add ray tracing to the gravity simulation to simulate gravitational lensing -2. Accretion disk : simulate accreciate disk using the ray tracing + the halos -3. Spacetime curvature : demonstrate visually the "trapdoor in spacetime" that is black holes using spacetime grid -4. [optional] try to make it run realtime ;D - -I hope it works :/ - - - -Edit: After completion of project - - -Thank you everyone for checking out the video, if you haven't it explains code in detail: https://www.youtube.com/watch?v=8-B6ryuBkCM - -Quickly to run the file you should install opengl (GLFW) and OpenGL wrangler GLEW: https://glew.sourceforge.net/ -I'm on windows and used msys2 for installation though. - -How the code works: - -for 2D: simple, just run 2D_lensing.cpp with the nessesary dependencies installed. - -for 3D: black_hole.cpp and geodesic.comp work together to run the simuation faster using GPU, essentially it sends over a UBO and geodesic.comp runs heavy calculations using that data. -should work with nessesary dependencies installed, however I have only run it on windows with my GPU so am not sure! - -LMK if you would like an in-depth explanation of how the code works aswell :) +Black hole simulation project: diff --git a/black_hole.cpp b/black_hole.cpp deleted file mode 100644 index 8999365..0000000 --- a/black_hole.cpp +++ /dev/null @@ -1,710 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#define _USE_MATH_DEFINES -#include -#include -#include -#include -#include -#include -#include -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -using namespace glm; -using namespace std; -using Clock = std::chrono::high_resolution_clock; - -// VARS -double lastPrintTime = 0.0; -int framesCount = 0; -double c = 299792458.0; -double G = 6.67430e-11; -struct Ray; -bool Gravity = false; - -struct Camera { - // Center the camera orbit on the black hole at (0, 0, 0) - vec3 target = vec3(0.0f, 0.0f, 0.0f); // Always look at the black hole center - float radius = 6.34194e10f; - float minRadius = 1e10f, maxRadius = 1e12f; - - float azimuth = 0.0f; - float elevation = M_PI / 2.0f; - - float orbitSpeed = 0.01f; - float panSpeed = 0.01f; - double zoomSpeed = 25e9f; - - bool dragging = false; - bool panning = false; - bool moving = false; // For compute shader optimization - double lastX = 0.0, lastY = 0.0; - - // Calculate camera position in world space - vec3 position() const { - float clampedElevation = clamp(elevation, 0.01f, float(M_PI) - 0.01f); - // Orbit around (0,0,0) always - return vec3( - radius * sin(clampedElevation) * cos(azimuth), - radius * cos(clampedElevation), - radius * sin(clampedElevation) * sin(azimuth) - ); - } - void update() { - // Always keep target at black hole center - target = vec3(0.0f, 0.0f, 0.0f); - if(dragging | panning) { - moving = true; - } else { - moving = false; - } - } - - void processMouseMove(double x, double y) { - float dx = float(x - lastX); - float dy = float(y - lastY); - - if (dragging && panning) { - // Pan: Shift + Left or Middle Mouse - // Disable panning to keep camera centered on black hole - } - else if (dragging && !panning) { - // Orbit: Left mouse only - azimuth += dx * orbitSpeed; - elevation -= dy * orbitSpeed; - elevation = clamp(elevation, 0.01f, float(M_PI) - 0.01f); - } - - lastX = x; - lastY = y; - update(); - } - void processMouseButton(int button, int action, int mods, GLFWwindow* win) { - if (button == GLFW_MOUSE_BUTTON_LEFT || button == GLFW_MOUSE_BUTTON_MIDDLE) { - if (action == GLFW_PRESS) { - dragging = true; - // Disable panning so camera always orbits center - panning = false; - glfwGetCursorPos(win, &lastX, &lastY); - } else if (action == GLFW_RELEASE) { - dragging = false; - panning = false; - } - } - if (button == GLFW_MOUSE_BUTTON_RIGHT) { - if (action == GLFW_PRESS) { - Gravity = true; - } else if (action == GLFW_RELEASE) { - Gravity = false; - } - } - } - void processScroll(double xoffset, double yoffset) { - radius -= yoffset * zoomSpeed; - radius = clamp(radius, minRadius, maxRadius); - update(); - } - void processKey(int key, int scancode, int action, int mods) { - if (action == GLFW_PRESS && key == GLFW_KEY_G) { - Gravity = !Gravity; - cout << "[INFO] Gravity turned " << (Gravity ? "ON" : "OFF") << endl; - } - } -}; -Camera camera; - -struct BlackHole { - vec3 position; - double mass; - double radius; - double r_s; - - BlackHole(vec3 pos, float m) : position(pos), mass(m) {r_s = 2.0 * G * mass / (c*c);} - bool Intercept(float px, float py, float pz) const { - double dx = double(px) - double(position.x); - double dy = double(py) - double(position.y); - double dz = double(pz) - double(position.z); - double dist2 = dx * dx + dy * dy + dz * dz; - return dist2 < r_s * r_s; - } -}; -BlackHole SagA(vec3(0.0f, 0.0f, 0.0f), 8.54e36); // Sagittarius A black hole -struct ObjectData { - vec4 posRadius; // xyz = position, w = radius - vec4 color; // rgb = color, a = unused - float mass; - vec3 velocity = vec3(0.0f, 0.0f, 0.0f); // Initial velocity -}; -vector objects = { - { vec4(4e11f, 0.0f, 0.0f, 4e10f) , vec4(1,1,0,1), 1.98892e30 }, - { vec4(0.0f, 0.0f, 4e11f, 4e10f) , vec4(1,0,0,1), 1.98892e30 }, - { vec4(0.0f, 0.0f, 0.0f, SagA.r_s) , vec4(0,0,0,1), SagA.mass }, - //{ vec4(6e10f, 0.0f, 0.0f, 5e10f), vec4(0,1,0,1) } -}; - -struct Engine { - GLuint gridShaderProgram; - // -- Quad & Texture render -- // - GLFWwindow* window; - GLuint quadVAO; - GLuint texture; - GLuint shaderProgram; - GLuint computeProgram = 0; - // -- UBOs -- // - GLuint cameraUBO = 0; - GLuint diskUBO = 0; - GLuint objectsUBO = 0; - // -- grid mess vars -- // - GLuint gridVAO = 0; - GLuint gridVBO = 0; - GLuint gridEBO = 0; - int gridIndexCount = 0; - - int WIDTH = 800; // Window width - int HEIGHT = 600; // Window height - int COMPUTE_WIDTH = 200; // Compute resolution width - int COMPUTE_HEIGHT = 150; // Compute resolution height - float width = 100000000000.0f; // Width of the viewport in meters - float height = 75000000000.0f; // Height of the viewport in meters - - Engine() { - if (!glfwInit()) { - cerr << "GLFW init failed\n"; - exit(EXIT_FAILURE); - } - glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4); - glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3); - glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE); - window = glfwCreateWindow(WIDTH, HEIGHT, "Black Hole", nullptr, nullptr); - if (!window) { - cerr << "Failed to create GLFW window\n"; - glfwTerminate(); - exit(EXIT_FAILURE); - } - glfwMakeContextCurrent(window); - glewExperimental = GL_TRUE; - GLenum glewErr = glewInit(); - if (glewErr != GLEW_OK) { - cerr << "Failed to initialize GLEW: " - << (const char*)glewGetErrorString(glewErr) - << "\n"; - glfwTerminate(); - exit(EXIT_FAILURE); - } - cout << "OpenGL " << glGetString(GL_VERSION) << "\n"; - this->shaderProgram = CreateShaderProgram(); - gridShaderProgram = CreateShaderProgram("grid.vert", "grid.frag"); - - computeProgram = CreateComputeProgram("geodesic.comp"); - glGenBuffers(1, &cameraUBO); - glBindBuffer(GL_UNIFORM_BUFFER, cameraUBO); - glBufferData(GL_UNIFORM_BUFFER, 128, nullptr, GL_DYNAMIC_DRAW); // alloc ~128 bytes - glBindBufferBase(GL_UNIFORM_BUFFER, 1, cameraUBO); // binding = 1 matches shader - - glGenBuffers(1, &diskUBO); - glBindBuffer(GL_UNIFORM_BUFFER, diskUBO); - glBufferData(GL_UNIFORM_BUFFER, sizeof(float) * 4, nullptr, GL_DYNAMIC_DRAW); // 3 values + 1 padding - glBindBufferBase(GL_UNIFORM_BUFFER, 2, diskUBO); // binding = 2 matches compute shader - - glGenBuffers(1, &objectsUBO); - glBindBuffer(GL_UNIFORM_BUFFER, objectsUBO); - // allocate space for 16 objects: - // sizeof(int) + padding + 16×(vec4 posRadius + vec4 color) - GLsizeiptr objUBOSize = sizeof(int) + 3 * sizeof(float) - + 16 * (sizeof(vec4) + sizeof(vec4)) - + 16 * sizeof(float); // 16 floats for mass - glBufferData(GL_UNIFORM_BUFFER, objUBOSize, nullptr, GL_DYNAMIC_DRAW); - glBindBufferBase(GL_UNIFORM_BUFFER, 3, objectsUBO); // binding = 3 matches shader - - auto result = QuadVAO(); - this->quadVAO = result[0]; - this->texture = result[1]; - } - void generateGrid(const vector& objects) { - const int gridSize = 25; - const float spacing = 1e10f; // tweak this - - vector vertices; - vector indices; - - for (int z = 0; z <= gridSize; ++z) { - for (int x = 0; x <= gridSize; ++x) { - float worldX = (x - gridSize / 2) * spacing; - float worldZ = (z - gridSize / 2) * spacing; - - float y = 0.0f; - - // ✅ Warp grid using Schwarzschild geometry - for (const auto& obj : objects) { - vec3 objPos = vec3(obj.posRadius); - double mass = obj.mass; - double radius = obj.posRadius.w; - - double r_s = 2.0 * G * mass / (c * c); - double dx = worldX - objPos.x; - double dz = worldZ - objPos.z; - double dist = sqrt(dx * dx + dz * dz); - - // prevent sqrt of negative or divide-by-zero (inside or at the black hole center) - if (dist > r_s) { - double deltaY = 2.0 * sqrt(r_s * (dist - r_s)); - y += static_cast(deltaY) - 3e10f; - } else { - // 🔴 For points inside or at r_s: make it dip down sharply - y += 2.0f * static_cast(sqrt(r_s * r_s)) - 3e10f; // or add a deep pit - } - } - - vertices.emplace_back(worldX, y, worldZ); - } - } - - // 🧩 Add indices for GL_LINE rendering - for (int z = 0; z < gridSize; ++z) { - for (int x = 0; x < gridSize; ++x) { - int i = z * (gridSize + 1) + x; - indices.push_back(i); - indices.push_back(i + 1); - - indices.push_back(i); - indices.push_back(i + gridSize + 1); - } - } - - // 🔌 Upload to GPU - if (gridVAO == 0) glGenVertexArrays(1, &gridVAO); - if (gridVBO == 0) glGenBuffers(1, &gridVBO); - if (gridEBO == 0) glGenBuffers(1, &gridEBO); - - glBindVertexArray(gridVAO); - - glBindBuffer(GL_ARRAY_BUFFER, gridVBO); - glBufferData(GL_ARRAY_BUFFER, vertices.size() * sizeof(vec3), vertices.data(), GL_DYNAMIC_DRAW); - - glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, gridEBO); - glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices.size() * sizeof(GLuint), indices.data(), GL_STATIC_DRAW); - - glEnableVertexAttribArray(0); // location = 0 - glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(vec3), (void*)0); - - gridIndexCount = indices.size(); - - glBindVertexArray(0); - } - void drawGrid(const mat4& viewProj) { - glUseProgram(gridShaderProgram); - glUniformMatrix4fv(glGetUniformLocation(gridShaderProgram, "viewProj"), - 1, GL_FALSE, glm::value_ptr(viewProj)); - glBindVertexArray(gridVAO); - - glDisable(GL_DEPTH_TEST); - glEnable(GL_BLEND); - glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); - - glDrawElements(GL_LINES, gridIndexCount, GL_UNSIGNED_INT, 0); - - glBindVertexArray(0); - glEnable(GL_DEPTH_TEST); - } - void drawFullScreenQuad() { - glUseProgram(shaderProgram); // fragment + vertex shader - glBindVertexArray(quadVAO); - - glActiveTexture(GL_TEXTURE0); - glBindTexture(GL_TEXTURE_2D, texture); - glUniform1i(glGetUniformLocation(shaderProgram, "screenTexture"), 0); - - glDisable(GL_DEPTH_TEST); // draw as background - glDrawArrays(GL_TRIANGLE_STRIP, 0, 6); // 2 triangles - glEnable(GL_DEPTH_TEST); - } - GLuint CreateShaderProgram(){ - const char* vertexShaderSource = R"( - #version 330 core - layout (location = 0) in vec2 aPos; // Changed to vec2 - layout (location = 1) in vec2 aTexCoord; - out vec2 TexCoord; - void main() { - gl_Position = vec4(aPos, 0.0, 1.0); // Explicit z=0 - TexCoord = aTexCoord; - })"; - - const char* fragmentShaderSource = R"( - #version 330 core - in vec2 TexCoord; - out vec4 FragColor; - uniform sampler2D screenTexture; - void main() { - FragColor = texture(screenTexture, TexCoord); - })"; - - // vertex shader - GLuint vertexShader = glCreateShader(GL_VERTEX_SHADER); - glShaderSource(vertexShader, 1, &vertexShaderSource, nullptr); - glCompileShader(vertexShader); - - // fragment shader - GLuint fragmentShader = glCreateShader(GL_FRAGMENT_SHADER); - glShaderSource(fragmentShader, 1, &fragmentShaderSource, nullptr); - glCompileShader(fragmentShader); - - GLuint shaderProgram = glCreateProgram(); - glAttachShader(shaderProgram, vertexShader); - glAttachShader(shaderProgram, fragmentShader); - glLinkProgram(shaderProgram); - - glDeleteShader(vertexShader); - glDeleteShader(fragmentShader); - - return shaderProgram; - }; - GLuint CreateShaderProgram(const char* vertPath, const char* fragPath) { - auto loadShader = [](const char* path, GLenum type) -> GLuint { - std::ifstream in(path); - if (!in.is_open()) { - std::cerr << "Failed to open shader: " << path << "\n"; - exit(EXIT_FAILURE); - } - std::stringstream ss; - ss << in.rdbuf(); - std::string srcStr = ss.str(); - const char* src = srcStr.c_str(); - - GLuint shader = glCreateShader(type); - glShaderSource(shader, 1, &src, nullptr); - glCompileShader(shader); - - GLint success; - glGetShaderiv(shader, GL_COMPILE_STATUS, &success); - if (!success) { - GLint logLen; - glGetShaderiv(shader, GL_INFO_LOG_LENGTH, &logLen); - std::vector log(logLen); - glGetShaderInfoLog(shader, logLen, nullptr, log.data()); - std::cerr << "Shader compile error (" << path << "):\n" << log.data() << "\n"; - exit(EXIT_FAILURE); - } - return shader; - }; - - GLuint vertShader = loadShader(vertPath, GL_VERTEX_SHADER); - GLuint fragShader = loadShader(fragPath, GL_FRAGMENT_SHADER); - - GLuint program = glCreateProgram(); - glAttachShader(program, vertShader); - glAttachShader(program, fragShader); - glLinkProgram(program); - - GLint linkSuccess; - glGetProgramiv(program, GL_LINK_STATUS, &linkSuccess); - if (!linkSuccess) { - GLint logLen; - glGetProgramiv(program, GL_INFO_LOG_LENGTH, &logLen); - std::vector log(logLen); - glGetProgramInfoLog(program, logLen, nullptr, log.data()); - std::cerr << "Shader link error:\n" << log.data() << "\n"; - exit(EXIT_FAILURE); - } - - glDeleteShader(vertShader); - glDeleteShader(fragShader); - - return program; - } - GLuint CreateComputeProgram(const char* path) { - // 1) read GLSL source - std::ifstream in(path); - if(!in.is_open()) { - std::cerr << "Failed to open compute shader: " << path << "\n"; - exit(EXIT_FAILURE); - } - std::stringstream ss; - ss << in.rdbuf(); - std::string srcStr = ss.str(); - const char* src = srcStr.c_str(); - - // 2) compile - GLuint cs = glCreateShader(GL_COMPUTE_SHADER); - glShaderSource(cs, 1, &src, nullptr); - glCompileShader(cs); - GLint ok; - glGetShaderiv(cs, GL_COMPILE_STATUS, &ok); - if(!ok) { - GLint logLen; - glGetShaderiv(cs, GL_INFO_LOG_LENGTH, &logLen); - std::vector log(logLen); - glGetShaderInfoLog(cs, logLen, nullptr, log.data()); - std::cerr << "Compute shader compile error:\n" << log.data() << "\n"; - exit(EXIT_FAILURE); - } - - // 3) link - GLuint prog = glCreateProgram(); - glAttachShader(prog, cs); - glLinkProgram(prog); - glGetProgramiv(prog, GL_LINK_STATUS, &ok); - if(!ok) { - GLint logLen; - glGetProgramiv(prog, GL_INFO_LOG_LENGTH, &logLen); - std::vector log(logLen); - glGetProgramInfoLog(prog, logLen, nullptr, log.data()); - std::cerr << "Compute shader link error:\n" << log.data() << "\n"; - exit(EXIT_FAILURE); - } - - glDeleteShader(cs); - return prog; - } - void dispatchCompute(const Camera& cam) { - // determine target compute‐res - int cw = cam.moving ? COMPUTE_WIDTH : 200; - int ch = cam.moving ? COMPUTE_HEIGHT : 150; - - // 1) reallocate the texture if needed - glBindTexture(GL_TEXTURE_2D, texture); - glTexImage2D(GL_TEXTURE_2D, - 0, // mip - GL_RGBA8, // internal format - cw, // width - ch, // height - 0, GL_RGBA, - GL_UNSIGNED_BYTE, - nullptr); - - // 2) bind compute program & UBOs - glUseProgram(computeProgram); - uploadCameraUBO(cam); - uploadDiskUBO(); - uploadObjectsUBO(objects); - - // 3) bind it as image unit 0 - glBindImageTexture(0, texture, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA8); - - // 4) dispatch grid - GLuint groupsX = (GLuint)std::ceil(cw / 16.0f); - GLuint groupsY = (GLuint)std::ceil(ch / 16.0f); - glDispatchCompute(groupsX, groupsY, 1); - - // 5) sync - glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT); - } - void uploadCameraUBO(const Camera& cam) { - struct UBOData { - vec3 pos; float _pad0; - vec3 right; float _pad1; - vec3 up; float _pad2; - vec3 forward; float _pad3; - float tanHalfFov; - float aspect; - bool moving; - int _pad4; - } data; - vec3 fwd = normalize(cam.target - cam.position()); - vec3 up = vec3(0, 1, 0); // y axis is up, so disk is in x-z plane - vec3 right = normalize(cross(fwd, up)); - up = cross(right, fwd); - - data.pos = cam.position(); - data.right = right; - data.up = up; - data.forward = fwd; - data.tanHalfFov = tan(radians(60.0f * 0.5f)); - data.aspect = float(WIDTH) / float(HEIGHT); - data.moving = cam.dragging || cam.panning; - - glBindBuffer(GL_UNIFORM_BUFFER, cameraUBO); - glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(UBOData), &data); - } - void uploadObjectsUBO(const vector& objs) { - struct UBOData { - int numObjects; - float _pad0, _pad1, _pad2; // <-- pad out to 16 bytes - vec4 posRadius[16]; - vec4 color[16]; - float mass[16]; - } data; - - size_t count = std::min(objs.size(), size_t(16)); - data.numObjects = static_cast(count); - - for (size_t i = 0; i < count; ++i) { - data.posRadius[i] = objs[i].posRadius; - data.color[i] = objs[i].color; - data.mass[i] = objs[i].mass; - } - - // Upload - glBindBuffer(GL_UNIFORM_BUFFER, objectsUBO); - glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(data), &data); - } - void uploadDiskUBO() { - // disk - float r1 = SagA.r_s * 2.2f; // inner radius just outside the event horizon - float r2 = SagA.r_s * 5.2f; // outer radius of the disk - float num = 2.0; // number of rays - float thickness = 1e9f; // padding for std140 alignment - float diskData[4] = { r1, r2, num, thickness }; - - glBindBuffer(GL_UNIFORM_BUFFER, diskUBO); - glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(diskData), diskData); - } - - vector QuadVAO(){ - float quadVertices[] = { - // positions // texCoords - -1.0f, 1.0f, 0.0f, 1.0f, // top left - -1.0f, -1.0f, 0.0f, 0.0f, // bottom left - 1.0f, -1.0f, 1.0f, 0.0f, // bottom right - - -1.0f, 1.0f, 0.0f, 1.0f, // top left - 1.0f, -1.0f, 1.0f, 0.0f, // bottom right - 1.0f, 1.0f, 1.0f, 1.0f // top right - }; - - GLuint VAO, VBO; - glGenVertexArrays(1, &VAO); - glGenBuffers(1, &VBO); - - glBindVertexArray(VAO); - glBindBuffer(GL_ARRAY_BUFFER, VBO); - glBufferData(GL_ARRAY_BUFFER, sizeof(quadVertices), quadVertices, GL_STATIC_DRAW); - - glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)0); - glEnableVertexAttribArray(0); - glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)(2 * sizeof(float))); - glEnableVertexAttribArray(1); - - GLuint texture; - glGenTextures(1, &texture); - glBindTexture(GL_TEXTURE_2D, texture); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - glBindTexture(GL_TEXTURE_2D, texture); - glTexImage2D(GL_TEXTURE_2D, - 0, // mip - GL_RGBA8, // internal format - COMPUTE_WIDTH, - COMPUTE_HEIGHT, - 0, - GL_RGBA, - GL_UNSIGNED_BYTE, - nullptr); - vector VAOtexture = {VAO, texture}; - return VAOtexture; - } - void renderScene() { - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - glUseProgram(shaderProgram); - glBindVertexArray(quadVAO); - // make sure your fragment shader samples from texture unit 0: - glActiveTexture(GL_TEXTURE0); - glBindTexture(GL_TEXTURE_2D, texture); - glDrawArrays(GL_TRIANGLES, 0, 6); - glfwSwapBuffers(window); - glfwPollEvents(); - }; -}; -Engine engine; -void setupCameraCallbacks(GLFWwindow* window) { - glfwSetWindowUserPointer(window, &camera); - - glfwSetMouseButtonCallback(window, [](GLFWwindow* win, int button, int action, int mods) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(win); - cam->processMouseButton(button, action, mods, win); - }); - - glfwSetCursorPosCallback(window, [](GLFWwindow* win, double x, double y) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(win); - cam->processMouseMove(x, y); - }); - - glfwSetScrollCallback(window, [](GLFWwindow* win, double xoffset, double yoffset) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(win); - cam->processScroll(xoffset, yoffset); - }); - - glfwSetKeyCallback(window, [](GLFWwindow* win, int key, int scancode, int action, int mods) { - Camera* cam = (Camera*)glfwGetWindowUserPointer(win); - cam->processKey(key, scancode, action, mods); - }); -} - - -// -- MAIN -- // -int main() { - setupCameraCallbacks(engine.window); - vector pixels(engine.WIDTH * engine.HEIGHT * 3); - - auto t0 = Clock::now(); - lastPrintTime = chrono::duration(t0.time_since_epoch()).count(); - - double lastTime = glfwGetTime(); - int renderW = 800, renderH = 600, numSteps = 80000; - while (!glfwWindowShouldClose(engine.window)) { - glClearColor(0.0f, 0.0f, 0.0f, 1.0f); // optional, but good practice - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - - double now = glfwGetTime(); - double dt = now - lastTime; // seconds since last frame - lastTime = now; - - // Gravity - for (auto& obj : objects) { - for (auto& obj2 : objects) { - if (&obj == &obj2) continue; // skip self-interaction - float dx = obj2.posRadius.x - obj.posRadius.x; - float dy = obj2.posRadius.y - obj.posRadius.y; - float dz = obj2.posRadius.z - obj.posRadius.z; - float distance = sqrt(dx * dx + dy * dy + dz * dz); - if (distance > 0) { - vector direction = {dx / distance, dy / distance, dz / distance}; - //distance *= 1000; - double Gforce = (G * obj.mass * obj2.mass) / (distance * distance); - - double acc1 = Gforce / obj.mass; - std::vector acc = {direction[0] * acc1, direction[1] * acc1, direction[2] * acc1}; - if (Gravity) { - obj.velocity.x += acc[0]; - obj.velocity.y += acc[1]; - obj.velocity.z += acc[2]; - - obj.posRadius.x += obj.velocity.x; - obj.posRadius.y += obj.velocity.y; - obj.posRadius.z += obj.velocity.z; - cout << "velocity: " <= disk_r1 && r <= disk_r2); -} - -void main() { - 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 - 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); - Ray ray = initRay(cam.camPos, dir); - - vec4 color = vec4(0.0); - vec3 prevPos = vec3(ray.x, ray.y, ray.z); - float lambda = 0.0; - - bool hitBlackHole = false; - bool hitDisk = false; - bool hitObject = false; - - int steps = cam.moving ? 60000 : 60000; - - for (int i = 0; i < steps; ++i) { - 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)) { hitDisk = true; break; } - if (interceptObject(ray)) { hitObject = true; break; } - prevPos = newPos; - if (ray.r > ESCAPE_R) break; - } - - if (hitDisk) { - 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) { - // 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); - } - - imageStore(outImage, pix, color); -} diff --git a/ray_tracing.cpp b/ray_tracing.cpp index 4479c30..9633c3e 100644 --- a/ray_tracing.cpp +++ b/ray_tracing.cpp @@ -5,286 +5,58 @@ #include #include #include -#include using namespace glm; -// global vars -const int WIDTH = 800; -const int HEIGHT = 600; - -// functions - -// structures and classes :D -class Engine{ -public: - // -- Quad & Texture render - GLFWwindow* window; - GLuint quadVAO; - GLuint texture; - GLuint shaderProgram; - - Engine(){ - this->window = StartGLFW(); - this->shaderProgram = CreateShaderProgram(); - - auto result = QuadVAO(); - this->quadVAO = result[0]; - this->texture = result[1]; - } - GLFWwindow* StartGLFW(){ - if(!glfwInit()){ - std::cerr<<"glfw failed init, PANIC PANIC!"< QuadVAO(){ - float quadVertices[] = { - // positions // texCoords - -1.0f, 1.0f, 0.0f, 1.0f, // top left - -1.0f, -1.0f, 0.0f, 0.0f, // bottom left - 1.0f, -1.0f, 1.0f, 0.0f, // bottom right - - -1.0f, 1.0f, 0.0f, 1.0f, // top left - 1.0f, -1.0f, 1.0f, 0.0f, // bottom right - 1.0f, 1.0f, 1.0f, 1.0f // top right - - }; - - GLuint VAO, VBO; - glGenVertexArrays(1, &VAO); - glGenBuffers(1, &VBO); - - glBindVertexArray(VAO); - glBindBuffer(GL_ARRAY_BUFFER, VBO); - glBufferData(GL_ARRAY_BUFFER, sizeof(quadVertices), quadVertices, GL_STATIC_DRAW); - - glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)0); - glEnableVertexAttribArray(0); - glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)(2 * sizeof(float))); - glEnableVertexAttribArray(1); - - GLuint texture; - glGenTextures(1, &texture); - glBindTexture(GL_TEXTURE_2D, texture); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - std::vector VAOtexture = {VAO, texture}; - return VAOtexture; - } - void renderScene(std::vector pixels) { - // update texture w/ ray-tracing results - glBindTexture(GL_TEXTURE_2D, texture); - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, WIDTH, HEIGHT, 0, GL_RGB, - GL_UNSIGNED_BYTE, pixels.data()); - - // clear screen and draw textured quad - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); - glUseProgram(shaderProgram); - - GLint textureLocation = glGetUniformLocation(shaderProgram, "screenTexture"); - glUniform1i(textureLocation, 0); - - glBindVertexArray(quadVAO); - glDrawArrays(GL_TRIANGLES, 0, 6); - - glfwSwapBuffers(window); - glfwPollEvents(); - }; -}; - +GLFWwindow* StartGLU(); struct Ray{ vec3 direction; vec3 origin; - Ray(vec3 o, vec3 d) : origin(o), direction(normalize(d)){} -}; -struct Material{ - vec3 color; - float specular; - float emission; - Material(vec3 c, float s, float e) : color(c), specular(s), emission(e) {} + Ray(vec3 direction, vec3 origin) : direction(direction), origin(origin) {} }; struct Object{ - vec3 centre; float radius; - Material material; - - Object(vec3 c, float r, Material m) : centre(c), radius(r), material(m) {} - - bool Intersect(Ray &ray, float &t){ - vec3 oc = ray.origin - centre; - float a = glm::dot(ray.direction, ray.direction); // ray direction scale by t - float b = 2.0f * glm::dot(oc, ray.direction); // - float c = glm::dot(oc, oc) - radius * radius; // adjustment by sphere radius - double discriminant = b*b - 4*a*c; - if(discriminant < 0){return false;} // no intersection with sphere - - float intercept = (-b - sqrt(discriminant)) / (2.0f*a); - if(intercept < 0){ - intercept = (-b + sqrt(discriminant)) / (2.0f*a); - if(intercept<0){return false;} // intersection is behind origin - } - t = intercept; - return true; - }; - - vec3 getNormal(vec3 &point) const{ - return normalize(point - centre); - } + vec3 centre; + Object(float radius, vec3 centre) : centre(centre), radius(radius){} }; -class Scene { -public: - std::vector objs; - vec3 lightPos; - Scene() : lightPos(5.0f, 5.0f, 5.0f) {} +int main() { + GLFWwindow* window = StartGLU(); - vec3 trace(Ray &ray){ - float closest = INFINITY; - const Object* hitObj = nullptr; + while(!glfwWindowShouldClose(window)){ - for(auto& obj : objs){ - float t; // distance to intersection - if(obj.Intersect(ray, t)){ - if(t < closest) { - closest = t; - hitObj = &obj; - } - } - }; - if(hitObj){ - vec3 hitPoint = ray.origin + ray.direction * closest; // point on obj hit by ray - vec3 normal = hitObj->getNormal(hitPoint); - vec3 lightDir = normalize(lightPos - hitPoint); // direction light to hitpoint - - float diff = std::max(glm::dot(normal, lightDir), 0.0f); // diffuse lighting - - Ray shadowRay(hitPoint + normal * 0.001f, lightDir); // slightly up to avoid errors ;P - // check if is in shadow - bool inShadow = false; - - // Actually check for shadows by testing if any object blocks light - for(auto& obj : objs) { - float t; - if(obj.Intersect(shadowRay, t)) { - inShadow = true; - break; - } - } - - vec3 color = hitObj->material.color; - float ambient = 0.1f; // minimum light level - - if (inShadow) { - return color * ambient; - } - - return color * (ambient + diff * 0.9f); - } - - return vec3(0.0f, 0.0f, 0.1f); + glfwSwapBuffers(window); + glfwPollEvents(); } -}; - - -// --- main loop ---- // -int main(){ - Engine engine; - Scene scene; - - scene.objs = { - Object(vec3(0.0f, 0.0f, -5.0f), 2.0f, Material(vec3(1.0f, 0.2f, 0.2f), 0.5f, 0.0f)), // Moved further back and made bigger - Object(vec3(3.0f, 0.0f, -7.0f), 1.5f, Material(vec3(0.2f, 1.0f, 0.2f), 0.5f, 0.0f)) // Adjusted position and size - }; - // -- loop -- // - std::vector pixels(WIDTH * HEIGHT * 3); - while(!glfwWindowShouldClose(engine.window)){ - glClear(GL_COLOR_BUFFER_BIT); - - // render texture (pxl by pxl) - for(int y = 0; y < HEIGHT; ++y){ - for(int x = 0; x < WIDTH; ++x){ - float aspectRatio = float(WIDTH) / float(HEIGHT); - float u = float(x) / float(WIDTH); - float v = float(y) / float(HEIGHT); - - // direction of ray threw camera - vec3 direction( - (2.0f * u - 1.0f) * aspectRatio, - -(2.0f * v - 1.0f), // Flipped to correct orientation - -1.0f // Forward direction (negative z) - ); - Ray ray(vec3(0.0f, 0.0f, 0.0f), normalize(direction)); - vec3 color = scene.trace(ray); - - int index = (y * WIDTH + x) * 3; - pixels[index + 0] = static_cast(color.r * 255); - pixels[index + 1] = static_cast(color.g * 255); - pixels[index + 2] = static_cast(color.b * 255); - } - } - - engine.renderScene(pixels); - } - glfwTerminate(); } +GLFWwindow* StartGLU() { + if (!glfwInit()) { + std::cout << "Failed to initialize GLFW, panic" << std::endl; + return nullptr; + } + GLFWwindow* window = glfwCreateWindow(800, 600, "RAY_TRACING", NULL, NULL); + if (!window) { + std::cerr << "Failed to create GLFW window, PANIC PANIC PANIC!" << std::endl; + glfwTerminate(); + return nullptr; + } + glfwMakeContextCurrent(window); -// func dec's + glewExperimental = GL_TRUE; + if (glewInit() != GLEW_OK) { + std::cerr << "Failed to initialize GLEW." << std::endl; + glfwTerminate(); + return nullptr; + } + + glEnable(GL_DEPTH_TEST); + glViewport(0, 0, 800, 600); + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // Standard blending for transparency + + return window; +} diff --git a/vs_code/c_cpp_properties.json b/vs_code/c_cpp_properties.json deleted file mode 100644 index 13ff85a..0000000 --- a/vs_code/c_cpp_properties.json +++ /dev/null @@ -1,18 +0,0 @@ -{ - "configurations": [ - { - "name": "Win32", - "includePath": [ - "${workspaceFolder}/**", - "C:/msys64/mingw64/include", - "C:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v12.8/include" - ], - "defines": [], - "compilerPath": "C:/msys64/mingw64/bin/g++.exe", - "cStandard": "c17", - "cppStandard": "c++17", - "intelliSenseMode": "windows-gcc-x64" - } - ], - "version": 4 -} \ No newline at end of file diff --git a/vs_code/launch.json b/vs_code/launch.json deleted file mode 100644 index 4d1a665..0000000 --- a/vs_code/launch.json +++ /dev/null @@ -1,47 +0,0 @@ -{ - "version": "0.2.0", - "configurations": [ - { - "name": "g++ Debug", - "type": "cppdbg", - "request": "launch", - "program": "${workspaceFolder}/bin/black_hole.exe", - "args": [], - "stopAtEntry": false, - "cwd": "${workspaceFolder}/src", - "environment": [], - "externalConsole": false, - "MIMode": "gdb", - "miDebuggerPath": "C:/msys64/mingw64/bin/gdb.exe", - "setupCommands": [ - { - "description": "Enable pretty-printing for gdb", - "text": "-enable-pretty-printing", - "ignoreFailures": true - } - ], - "preLaunchTask": "build" - }, - { - "name": "CUDA Debug", - "type": "cppdbg", - "request": "launch", - "program": "${workspaceFolder}/bin/black_hole.exe", - "args": [], - "stopAtEntry": false, - "cwd": "${workspaceFolder}/src", - "environment": [], - "externalConsole": false, - "MIMode": "gdb", - "miDebuggerPath": "C:/msys64/mingw64/bin/gdb.exe", - "setupCommands": [ - { - "description": "Enable pretty-printing for gdb", - "text": "-enable-pretty-printing", - "ignoreFailures": true - } - ], - "preLaunchTask": "build-cuda" - } - ] -} \ No newline at end of file diff --git a/vs_code/settings.json b/vs_code/settings.json deleted file mode 100644 index 145c23d..0000000 --- a/vs_code/settings.json +++ /dev/null @@ -1,77 +0,0 @@ -{ - "files.associations": { - "ostream": "cpp", - "array": "cpp", - "atomic": "cpp", - "bit": "cpp", - "*.tcc": "cpp", - "cctype": "cpp", - "charconv": "cpp", - "clocale": "cpp", - "cmath": "cpp", - "compare": "cpp", - "concepts": "cpp", - "cstdarg": "cpp", - "cstddef": "cpp", - "cstdint": "cpp", - "cstdio": "cpp", - "cstdlib": "cpp", - "ctime": "cpp", - "cwchar": "cpp", - "cwctype": "cpp", - "deque": "cpp", - "string": "cpp", - "unordered_map": "cpp", - "vector": "cpp", - "exception": "cpp", - "algorithm": "cpp", - "functional": "cpp", - "iterator": "cpp", - "memory": "cpp", - "memory_resource": "cpp", - "numeric": "cpp", - "optional": "cpp", - "random": "cpp", - "string_view": "cpp", - "system_error": "cpp", - "tuple": "cpp", - "type_traits": "cpp", - "utility": "cpp", - "format": "cpp", - "initializer_list": "cpp", - "iosfwd": "cpp", - "iostream": "cpp", - "istream": "cpp", - "limits": "cpp", - "new": "cpp", - "numbers": "cpp", - "span": "cpp", - "stdexcept": "cpp", - "streambuf": "cpp", - "text_encoding": "cpp", - "cinttypes": "cpp", - "typeinfo": "cpp", - "variant": "cpp", - "chrono": "cpp", - "ratio": "cpp", - "iomanip": "cpp", - "semaphore": "cpp", - "sstream": "cpp", - "stop_token": "cpp", - "thread": "cpp", - "cstring": "cpp", - "ios": "cpp", - "list": "cpp", - "xfacet": "cpp", - "xhash": "cpp", - "xiosbase": "cpp", - "xlocale": "cpp", - "xlocinfo": "cpp", - "xlocnum": "cpp", - "xmemory": "cpp", - "xstddef": "cpp", - "xstring": "cpp", - "xtr1common": "cpp", - "xutility": "cpp" - } -} \ No newline at end of file diff --git a/vs_code/tasks.json b/vs_code/tasks.json deleted file mode 100644 index b9a6a36..0000000 --- a/vs_code/tasks.json +++ /dev/null @@ -1,57 +0,0 @@ -{ - "version": "2.0.0", - "tasks": [ - { - "label": "build", - "type": "cppbuild", - "command": "C:/msys64/mingw64/bin/g++.exe", - "args": [ - "-fdiagnostics-color=always", - "-g", - "${workspaceFolder}/src/black_hole.cpp", - "-o", - "${workspaceFolder}/bin/black_hole.exe", - "-IC:/msys64/mingw64/include", - "-LC:/msys64/mingw64/lib", - "-lglfw3", - "-lglew32", - "-lopengl32", - "-lgdi32" - ], - "options": { - "cwd": "${workspaceFolder}" - }, - "problemMatcher": [ - "$gcc" - ], - "group": { - "kind": "build", - "isDefault": true - }, - "detail": "compiler: C:/msys64/mingw64/bin/g++.exe" - }, - { - "label": "build-cuda", - "type": "shell", - "command": "\"C:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v12.8/bin/nvcc.exe\"", - "args": [ - "-o", - "${workspaceFolder}/bin/black_hole_cuda.exe", - "${workspaceFolder}/src/black_hole.cu", - "-IC:/msys64/mingw64/include", - "-IC:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v12.8/include", - "-LC:/msys64/mingw64/lib", - "-LC:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v12.8/lib/x64", - "-lglfw3", - "-lglew32", - "-lopengl32", - "-lgdi32", - "-lcudart" - ], - "options": { - "cwd": "${workspaceFolder}" - }, - "group": "build" - } - ] -} \ No newline at end of file