Completed Black Hole README.md

while moving lower res

2D_lensing.cpp

@@ -0,0 +1,231 @@
+#include <GL/glew.h>
+#include <GLFW/glfw3.h>
+#include <glm/glm.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+#include <glm/gtc/type_ptr.hpp>
+#include <vector>
+#include <iostream>
+#define _USE_MATH_DEFINES
+#include <cmath>
+#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<vec2> 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<Ray>& 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<Ray> 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;
+}

CPU-geodesic.cpp

@@ -0,0 +1,496 @@
+#include <GL/glew.h>
+#include <GLFW/glfw3.h>
+#include <glm/glm.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+#include <glm/gtc/type_ptr.hpp>
+#include <vector>
+#include <iostream>
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <sstream>
+#include <iomanip>
+#include <cstring>
+#include <chrono>
+#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<GLuint> 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<GLuint> VAOtexture = {VAO, texture};
+        return VAOtexture;
+    }
+    void renderScene(const vector<unsigned char>& 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<unsigned char>& 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<unsigned char> pixels(engine.WIDTH * engine.HEIGHT * 3);
+
+    auto t0 = Clock::now();
+    lastPrintTime = std::chrono::duration<double>(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<double>(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<double>(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);

README.md

@@ -1,2 +1,29 @@
 # black_hole
-Black hole simulation project:
+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.cpp

@@ -0,0 +1,710 @@
+#include <GL/glew.h>
+#include <GLFW/glfw3.h>
+#include <glm/glm.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+#include <glm/gtc/type_ptr.hpp>
+#include <vector>
+#include <iostream>
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <sstream>
+#include <iomanip>
+#include <cstring>
+#include <chrono>
+#include <fstream>
+#include <sstream>
+#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<ObjectData> 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<ObjectData>& objects) {
+        const int gridSize = 25;
+        const float spacing = 1e10f;  // tweak this
+
+        vector<vec3> vertices;
+        vector<GLuint> 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<float>(deltaY) - 3e10f;
+                    } else {
+                        // 🔴 For points inside or at r_s: make it dip down sharply
+                        y += 2.0f * static_cast<float>(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<char> 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<char> 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<char> 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<char> 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<ObjectData>& 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<int>(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<GLuint> 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<GLuint> 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<unsigned char> pixels(engine.WIDTH * engine.HEIGHT * 3);
+
+    auto t0 = Clock::now();
+    lastPrintTime = chrono::duration<double>(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<double> 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<double> 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: " <<obj.velocity.x<<", " <<obj.velocity.y<<", " <<obj.velocity.z<<endl;
+                        }
+                    }
+            }
+        }
+
+
+
+        // ---------- GRID ------------- //
+        // 2) rebuild grid mesh on CPU
+        engine.generateGrid(objects);
+        // 5) overlay the bent grid
+        mat4 view = lookAt(camera.position(), camera.target, vec3(0,1,0));
+        mat4 proj = perspective(radians(60.0f), float(engine.COMPUTE_WIDTH)/engine.COMPUTE_HEIGHT, 1e9f, 1e14f);
+        mat4 viewProj = proj * view;
+        engine.drawGrid(viewProj);
+
+        // ---------- RUN RAYTRACER ------------- //
+        glViewport(0, 0, engine.WIDTH, engine.HEIGHT);
+        engine.dispatchCompute(camera);
+        engine.drawFullScreenQuad();
+
+        // 6) present to screen
+        glfwSwapBuffers(engine.window);
+        glfwPollEvents();
+    }
+
+    glfwDestroyWindow(engine.window);
+    glfwTerminate();
+    return 0;
+}

geodesic.comp

@@ -0,0 +1,177 @@
+#version 430
+layout(local_size_x = 16, local_size_y = 16) in;
+
+layout(binding = 0, rgba8) writeonly uniform image2D outImage;
+layout(std140, binding = 1) uniform Camera {
+    vec3 camPos;     float _pad0;
+    vec3 camRight;   float _pad1;
+    vec3 camUp;      float _pad2;
+    vec3 camForward; float _pad3;
+    float tanHalfFov;
+    float aspect;
+    bool moving;
+    int   _pad4;
+} cam;
+
+layout(std140, binding = 2) uniform Disk {
+    float disk_r1;
+    float disk_r2;
+    float disk_num;
+    float thickness;
+};
+
+layout(std140, binding = 3) uniform Objects {
+    int numObjects;
+    vec4 objPosRadius[16];
+    vec4 objColor[16];
+    float  mass[16]; 
+};
+
+const float SagA_rs = 1.269e10;
+const float D_LAMBDA = 1e7;
+const double ESCAPE_R = 1e30;
+
+// 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;
+};
+Ray initRay(vec3 pos, vec3 dir) {
+    Ray ray;
+    ray.x = pos.x; ray.y = pos.y; ray.z = pos.z;
+    ray.r = length(pos);
+    ray.theta = acos(pos.z / ray.r);
+    ray.phi = atan(pos.y, pos.x);
+
+    float dx = dir.x, dy = dir.y, dz = dir.z;
+    ray.dr     = sin(ray.theta)*cos(ray.phi)*dx + sin(ray.theta)*sin(ray.phi)*dy + cos(ray.theta)*dz;
+    ray.dtheta = (cos(ray.theta)*cos(ray.phi)*dx + cos(ray.theta)*sin(ray.phi)*dy - sin(ray.theta)*dz) / ray.r;
+    ray.dphi   = (-sin(ray.phi)*dx + cos(ray.phi)*dy) / (ray.r * sin(ray.theta));
+
+    ray.L = ray.r * ray.r * sin(ray.theta) * ray.dphi;
+    float f = 1.0 - SagA_rs / ray.r;
+    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;
+}
+
+bool intercept(Ray ray, float 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 center = objPosRadius[i].xyz;
+        float radius = objPosRadius[i].w;
+        if (distance(P, center) <= radius) {
+            objectColor = objColor[i];
+            hitCenter = center;
+            hitRadius = radius;
+            return true;
+        }
+    }
+    return false;
+}
+
+void geodesicRHS(Ray ray, out vec3 d1, out vec3 d2) {
+    float r = ray.r, theta = ray.theta;
+    float dr = ray.dr, dtheta = ray.dtheta, dphi = ray.dphi;
+    float f = 1.0 - SagA_rs / r;
+    float dt_dL = ray.E / f;
+
+    d1 = vec3(dr, dtheta, dphi);
+    d2.x = - (SagA_rs / (2.0 * r*r)) * f * dt_dL * dt_dL
+         + (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;
+}
+void rk4Step(inout Ray ray, float dL) {
+    vec3 k1a, k1b;
+    geodesicRHS(ray, k1a, k1b);
+
+    ray.r      += dL * k1a.x;
+    ray.theta  += dL * k1a.y;
+    ray.phi    += dL * k1a.z;
+    ray.dr     += dL * k1b.x;
+    ray.dtheta += dL * k1b.y;
+    ray.dphi   += dL * k1b.z;
+
+    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);
+}
+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);
+}
+
+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);
+}

ray_tracing.cpp

@@ -5,58 +5,286 @@
 #include <glm/gtc/type_ptr.hpp>
 #include <vector>
 #include <iostream>
+#include <cmath>
 using namespace glm;
 
-GLFWwindow* StartGLU();
-struct Ray{
-    vec3 direction;
-    vec3 origin;
-    Ray(vec3 direction, vec3 origin) : direction(direction), origin(origin) {}
-};
-struct Object{
-    float radius;
-    vec3 centre;
-    Object(float radius, vec3 centre) : centre(centre), radius(radius){}
-};
+// global vars
+const int WIDTH = 800;
+const int HEIGHT = 600;
 
-int main() {
-    GLFWwindow* window = StartGLU();
+// functions
 
-    while(!glfwWindowShouldClose(window)){
+// 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!"<<std::endl;
+            return nullptr;
+        }
+
+        GLFWwindow* window = glfwCreateWindow(WIDTH, HEIGHT, "ray tracer", NULL, NULL);
+        glfwMakeContextCurrent(window);
+        
+        glewExperimental = GL_TRUE;
+        if (glewInit() != GLEW_OK) {
+            std::cerr << "Failed to initialize GLEW." << std::endl;
+            glfwTerminate();
+            return nullptr;
+        }
+
+        glViewport(0, 0, WIDTH, HEIGHT);
+        return window;
+    };
+    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;
+    };
+
+    std::vector<GLuint> 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<GLuint> VAOtexture = {VAO, texture};
+        return VAOtexture;
+    }
+    void renderScene(std::vector<unsigned char> 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();
+    };
+};
+
+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) {}
+};
+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);
     }
+};
+
+class Scene {
+public:
+    std::vector<Object> objs;
+    vec3 lightPos;
+    Scene() : lightPos(5.0f, 5.0f, 5.0f) {}
+
+    vec3 trace(Ray &ray){
+        float closest = INFINITY;
+        const Object* hitObj = nullptr;
+
+        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); 
+    }
+};
+
+
+// --- 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<unsigned char> 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<unsigned char>(color.r * 255);
+                pixels[index + 1] = static_cast<unsigned char>(color.g * 255);
+                pixels[index + 2] = static_cast<unsigned char>(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);
 
-    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;
-}
+// func dec's
 
 
 

vs_code/c_cpp_properties.json

@@ -0,0 +1,18 @@
+{
+    "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
+}

vs_code/launch.json

@@ -0,0 +1,47 @@
+{
+    "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"
+        }
+    ]
+}

vs_code/settings.json

@@ -0,0 +1,77 @@
+{
+    "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"
+    }
+}

vs_code/tasks.json

@@ -0,0 +1,57 @@
+{
+    "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"
+        }
+    ]
+}