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

@@ -9,3 +9,21 @@ I'm writing this as I'm beginning this project (hopefully I complete it ;D) here
 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 :)

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);
+}

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"
+        }
+    ]
+}