diff --git a/2D_lensing.cpp b/2D_lensing.cpp new file mode 100644 index 0000000..f2f4e79 --- /dev/null +++ b/2D_lensing.cpp @@ -0,0 +1,231 @@ +#include +#include +#include +#include +#include +#include +#include +#define _USE_MATH_DEFINES +#include +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +using namespace glm; +using namespace std; + +double c = 299792458.0; +double G = 6.67430e-11; + +struct Ray; +void rk4Step(Ray& ray, double dλ, double rs); + +// --- Structs --- // +struct Engine { + GLFWwindow* window; + int WIDTH = 800; + int HEIGHT = 600; + float width = 100000000000.0f; // Width of the viewport in meters + float height = 75000000000.0f; // Height of the viewport in meters + + // Navigation state + float offsetX = 0.0f, offsetY = 0.0f; + float zoom = 1.0f; + bool middleMousePressed = false; + double lastMouseX = 0.0, lastMouseY = 0; + + Engine() { + if (!glfwInit()) { + cerr << "Failed to initialize GLFW" << endl; + exit(EXIT_FAILURE); + } + window = glfwCreateWindow(WIDTH, HEIGHT, "Black Hole Simulation", NULL, NULL); + if (!window) { + cerr << "Failed to create GLFW window" << endl; + glfwTerminate(); + exit(EXIT_FAILURE); + } + glfwMakeContextCurrent(window); + glewExperimental = GL_TRUE; + if (glewInit() != GLEW_OK) { + cerr << "Failed to initialize GLEW" << endl; + glfwDestroyWindow(window); + glfwTerminate(); + exit(EXIT_FAILURE); + } + glViewport(0, 0, WIDTH, HEIGHT);; + } + + void run() { + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + double left = -width + offsetX; + double right = width + offsetX; + double bottom = -height + offsetY; + double top = height + offsetY; + glOrtho(left, right, bottom, top, -1.0, 1.0); + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + } +}; +Engine engine; +struct BlackHole { + vec3 position; + double mass; + double radius; + double r_s; + + BlackHole(vec3 pos, float m) : position(pos), mass(m) {r_s = 2.0 * G * mass / (c*c);} + void draw() { + glBegin(GL_TRIANGLE_FAN); + glColor3f(1.0f, 0.0f, 0.0f); // Red color for the black hole + glVertex2f(0.0f, 0.0f); // Center + for(int i = 0; i <= 100; i++) { + float angle = 2.0f * M_PI * i / 100; + float x = r_s * cos(angle); // Radius of 0.1 + float y = r_s * sin(angle); + glVertex2f(x, y); + } + glEnd(); + } +}; +BlackHole SagA(vec3(0.0f, 0.0f, 0.0f), 8.54e36); // Sagittarius A black hole +struct Ray{ + // -- cartesian coords -- // + double x; double y; + // -- polar coords -- // + double r; double phi; + double dr; double dphi; + vector trail; // trail of points + double E, L; // conserved quantities + + Ray(vec2 pos, vec2 dir) : x(pos.x), y(pos.y), r(sqrt(pos.x * pos.x + pos.y * pos.y)), phi(atan2(pos.y, pos.x)), dr(dir.x), dphi(dir.y) { + // step 1) get polar coords (r, phi) : + this->r = sqrt(x*x + y*y); + this->phi = atan2(y, x); + // step 2) seed velocities : + dr = dir.x * cos(phi) + dir.y * sin(phi); // m/s + dphi = ( -dir.x * sin(phi) + dir.y * cos(phi) ) / r; + // step 3) store conserved quantities + L = r*r * dphi; + double f = 1.0 - SagA.r_s/r; + double dt_dλ = sqrt( (dr*dr)/(f*f) + (r*r*dphi*dphi)/f ); + E = f * dt_dλ; + // step 4) start trail : + trail.push_back({x, y}); + } + void draw(const std::vector& rays) { + // draw current ray positions as points + glPointSize(2.0f); + glColor3f(1.0f, 0.0f, 0.0f); + glBegin(GL_POINTS); + for (const auto& ray : rays) { + glVertex2f(ray.x, ray.y); + } + glEnd(); + + // turn on blending for the trails + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + glLineWidth(1.0f); + + // draw each trail with fading alpha + for (const auto& ray : rays) { + size_t N = ray.trail.size(); + if (N < 2) continue; + + glBegin(GL_LINE_STRIP); + for (size_t i = 0; i < N; ++i) { + // older points (i=0) get alpha≈0, newer get alpha≈1 + float alpha = float(i) / float(N - 1); + glColor4f(1.0f, 1.0f, 1.0f, std::max(alpha, 0.05f)); + glVertex2f(ray.trail[i].x, ray.trail[i].y); + } + glEnd(); + } + + glDisable(GL_BLEND); + } + void step(double dλ, double rs) { + // 1) integrate (r,φ,dr,dφ) + if(r <= rs) return; // stop if inside the event horizon + rk4Step(*this, dλ, rs); + + // 2) convert back to cartesian x,y + x = r * cos(phi); + y = r * sin(phi); + + // 3) record the trail + trail.push_back({ float(x), float(y) }); + } +}; +vector rays; + +void geodesicRHS(const Ray& ray, double rhs[4], double rs) { + double r = ray.r; + double dr = ray.dr; + double dphi = ray.dphi; + double E = ray.E; + + double f = 1.0 - rs/r; + + // dr/dλ = dr + rhs[0] = dr; + // dφ/dλ = dphi + rhs[1] = dphi; + + // d²r/dλ² from Schwarzschild null geodesic: + double dt_dλ = E / f; + rhs[2] = + - (rs/(2*r*r)) * f * (dt_dλ*dt_dλ) + + (rs/(2*r*r*f)) * (dr*dr) + + (r - rs) * (dphi*dphi); + + // d²φ/dλ² = -2*(dr * dphi) / r + rhs[3] = -2.0 * dr * dphi / r; +} +void addState(const double a[4], const double b[4], double factor, double out[4]) { + for (int i = 0; i < 4; i++) + out[i] = a[i] + b[i] * factor; +} +void rk4Step(Ray& ray, double dλ, double rs) { + double y0[4] = { ray.r, ray.phi, ray.dr, ray.dphi }; + double k1[4], k2[4], k3[4], k4[4], temp[4]; + + geodesicRHS(ray, k1, rs); + addState(y0, k1, dλ/2.0, temp); + Ray r2 = ray; r2.r=temp[0]; r2.phi=temp[1]; r2.dr=temp[2]; r2.dphi=temp[3]; + geodesicRHS(r2, k2, rs); + + addState(y0, k2, dλ/2.0, temp); + Ray r3 = ray; r3.r=temp[0]; r3.phi=temp[1]; r3.dr=temp[2]; r3.dphi=temp[3]; + geodesicRHS(r3, k3, rs); + + addState(y0, k3, dλ, temp); + Ray r4 = ray; r4.r=temp[0]; r4.phi=temp[1]; r4.dr=temp[2]; r4.dphi=temp[3]; + geodesicRHS(r4, k4, rs); + + ray.r += (dλ/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); + ray.phi += (dλ/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]); + ray.dr += (dλ/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]); + ray.dphi += (dλ/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]); +} + + +int main () { + //rays.push_back(Ray(vec2(-1e11, 3.27606302719999999e10), vec2(c, 0.0f))); + while(!glfwWindowShouldClose(engine.window)) { + engine.run(); + SagA.draw(); + + for (auto& ray : rays) { + ray.step(1.0f, SagA.r_s); + ray.draw(rays); + } + + glfwSwapBuffers(engine.window); + glfwPollEvents(); + } + + return 0; +} diff --git a/CPU-geodesic.cpp b/CPU-geodesic.cpp new file mode 100644 index 0000000..e1f7e30 --- /dev/null +++ b/CPU-geodesic.cpp @@ -0,0 +1,496 @@ +#include +#include +#include +#include +#include +#include +#include +#define _USE_MATH_DEFINES +#include +#include +#include +#include +#include +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +using namespace glm; +using namespace std; +using Clock = std::chrono::high_resolution_clock; + +// VARS +double lastPrintTime = 0.0; +int framesCount = 0; +double c = 299792458.0; +double G = 6.67430e-11; +bool useGeodesics = false; + +struct Camera { + vec3 pos; + vec3 target; + float fovY; + float azimuth, elevation, radius; + float minRadius = 1e12f, maxRadius = 1e20f; + bool dragging = false; + bool panning = false; + double lastX = 0, lastY = 0; + + // Adjustable speeds + float orbitSpeed = 0.008f; + float panSpeed = 0.001f; + float zoomSpeed = 1.08f; // closer to 1 = slower zoom + + Camera() : azimuth(0), elevation(M_PI / 2.0f), radius(6.34194e10), fovY(60.0f) { + target = vec3(0, 0, 0); + updateVectors(); + } + + void updateVectors() { + pos.x = target.x + radius * sin(elevation) * cos(azimuth); + pos.y = target.y + radius * cos(elevation); + pos.z = target.z + radius * sin(elevation) * sin(azimuth); + } + void processMouse(GLFWwindow* window, double xpos, double ypos) { + float dx = float(xpos - lastX), dy = float(ypos - lastY); + if (dragging && !panning) { + // Orbit + azimuth -= dx * orbitSpeed; + elevation -= dy * orbitSpeed; + elevation = clamp(elevation, 0.01f, float(M_PI)-0.01f); + } else if (panning) { + // Pan (move target in camera plane) + vec3 forward = normalize(target - pos); + vec3 right = normalize(cross(forward, vec3(0,1,0))); + vec3 up = cross(right, forward); + target += -right * dx * panSpeed * radius + up * dy * panSpeed * radius; + } + updateVectors(); + lastX = xpos; lastY = ypos; + } + void processScroll(double yoffset) { + // Zoom (dolly in/out) + if (yoffset < 0) + radius *= pow(zoomSpeed, -yoffset); + else + radius /= pow(zoomSpeed, yoffset); + radius = clamp(radius, minRadius, maxRadius); + updateVectors(); + } + static void mouseButtonCallback(GLFWwindow* window, int button, int action, int mods) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(window); + if (button == GLFW_MOUSE_BUTTON_LEFT) { + if (action == GLFW_PRESS) { + cam->dragging = true; + cam->panning = (mods & GLFW_MOD_SHIFT); + double x, y; glfwGetCursorPos(window, &x, &y); + cam->lastX = x; cam->lastY = y; + } else if (action == GLFW_RELEASE) { + cam->dragging = false; + cam->panning = false; + } + } + } + static void cursorPosCallback(GLFWwindow* window, double xpos, double ypos) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(window); + cam->processMouse(window, xpos, ypos); + } + static void scrollCallback(GLFWwindow* window, double xoffset, double yoffset) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(window); + cam->processScroll(yoffset); + } +}; +Camera camera; + +struct Ray; +void rk4Step(Ray& ray, double dλ, double rs); + +struct Engine { + // -- Quad & Texture render -- // + GLFWwindow* window; + GLuint quadVAO; + GLuint texture; + GLuint shaderProgram; + int WIDTH = 800; + int HEIGHT = 600; + float width = 100000000000.0f; // Width of the viewport in meters + float height = 75000000000.0f; // Height of the viewport in meters + + Engine() { + if (!glfwInit()) { + cerr << "GLFW init failed\n"; + exit(EXIT_FAILURE); + } + glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3); + glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3); + glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE); + window = glfwCreateWindow(WIDTH, HEIGHT, "Black Hole", nullptr, nullptr); + if (!window) { + cerr << "Failed to create GLFW window\n"; + glfwTerminate(); + exit(EXIT_FAILURE); + } + glfwMakeContextCurrent(window); + glewExperimental = GL_TRUE; + GLenum glewErr = glewInit(); + if (glewErr != GLEW_OK) { + cerr << "Failed to initialize GLEW: " + << (const char*)glewGetErrorString(glewErr) + << "\n"; + glfwTerminate(); + exit(EXIT_FAILURE); + } + cout << "OpenGL " << glGetString(GL_VERSION) << "\n"; + this->shaderProgram = CreateShaderProgram(); + + auto result = QuadVAO(); + this->quadVAO = result[0]; + this->texture = result[1]; + } + GLuint CreateShaderProgram(){ + const char* vertexShaderSource = R"( + #version 330 core + layout (location = 0) in vec2 aPos; // Changed to vec2 + layout (location = 1) in vec2 aTexCoord; + out vec2 TexCoord; + void main() { + gl_Position = vec4(aPos, 0.0, 1.0); // Explicit z=0 + TexCoord = aTexCoord; + })"; + + const char* fragmentShaderSource = R"( + #version 330 core + in vec2 TexCoord; + out vec4 FragColor; + uniform sampler2D screenTexture; + void main() { + FragColor = texture(screenTexture, TexCoord); + })"; + + // vertex shader + GLuint vertexShader = glCreateShader(GL_VERTEX_SHADER); + glShaderSource(vertexShader, 1, &vertexShaderSource, nullptr); + glCompileShader(vertexShader); + + // fragment shader + GLuint fragmentShader = glCreateShader(GL_FRAGMENT_SHADER); + glShaderSource(fragmentShader, 1, &fragmentShaderSource, nullptr); + glCompileShader(fragmentShader); + + GLuint shaderProgram = glCreateProgram(); + glAttachShader(shaderProgram, vertexShader); + glAttachShader(shaderProgram, fragmentShader); + glLinkProgram(shaderProgram); + + glDeleteShader(vertexShader); + glDeleteShader(fragmentShader); + + return shaderProgram; + }; + vector QuadVAO(){ + float quadVertices[] = { + // positions // texCoords + -1.0f, 1.0f, 0.0f, 1.0f, // top left + -1.0f, -1.0f, 0.0f, 0.0f, // bottom left + 1.0f, -1.0f, 1.0f, 0.0f, // bottom right + + -1.0f, 1.0f, 0.0f, 1.0f, // top left + 1.0f, -1.0f, 1.0f, 0.0f, // bottom right + 1.0f, 1.0f, 1.0f, 1.0f // top right + + }; + + GLuint VAO, VBO; + glGenVertexArrays(1, &VAO); + glGenBuffers(1, &VBO); + + glBindVertexArray(VAO); + glBindBuffer(GL_ARRAY_BUFFER, VBO); + glBufferData(GL_ARRAY_BUFFER, sizeof(quadVertices), quadVertices, GL_STATIC_DRAW); + + glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)0); + glEnableVertexAttribArray(0); + glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 4 * sizeof(float), (void*)(2 * sizeof(float))); + glEnableVertexAttribArray(1); + + GLuint texture; + glGenTextures(1, &texture); + glBindTexture(GL_TEXTURE_2D, texture); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); + vector VAOtexture = {VAO, texture}; + return VAOtexture; + } + void renderScene(const vector& pixels, int texWidth, int texHeight) { + // update texture w/ ray-tracing results + glBindTexture(GL_TEXTURE_2D, texture); + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, texWidth, texHeight, 0, GL_RGB, GL_UNSIGNED_BYTE, pixels.data()); + + // clear screen and draw textured quad + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glUseProgram(shaderProgram); + + GLint textureLocation = glGetUniformLocation(shaderProgram, "screenTexture"); + glUniform1i(textureLocation, 0); + + glBindVertexArray(quadVAO); + glDrawArrays(GL_TRIANGLES, 0, 6); + + glfwSwapBuffers(window); + glfwPollEvents(); + }; + static void keyCallback(GLFWwindow* window, int key, int scancode, int action, int mods) { + if (action == GLFW_PRESS) { + if (key == GLFW_KEY_G) { + useGeodesics = !useGeodesics; + cout << "Geodesics: " << (useGeodesics ? "ON\n" : "OFF\n"); + } + } + } +}; +Engine engine; +struct BlackHole { + vec3 position; + double mass; + double radius; + double r_s; + + BlackHole(vec3 pos, float m) : position(pos), mass(m) {r_s = 2.0 * G * mass / (c*c);} + bool Intercept(float px, float py, float pz) const { + float dx = px - position.x; + float dy = py - position.y; + float dz = pz - position.z; + float dist2 = dx * dx + dy * dy + dz * dz; + return dist2 < r_s * r_s; + } +}; +BlackHole SagA(vec3(0.0f, 0.0f, 0.0f), 8.54e36); // Sagittarius A black hole +struct Ray{ + // -- cartesian coords -- // + double x; double y; double z; + // -- polar coords -- // + double r; double phi; double theta; + double dr; double dphi; double dtheta; + double E, L; // conserved quantities + + Ray(vec3 pos, vec3 dir) : x(pos.x), y(pos.y), z(pos.z) { + // Step 1: get spherical coords (r, theta, phi) + r = sqrt(x*x + y*y + z*z); + theta = acos(z / r); + phi = atan2(y, x); + + // Step 2: seed velocities (dr, dtheta, dphi) + // Convert direction to spherical basis + double dx = dir.x, dy = dir.y, dz = dir.z; + dr = sin(theta)*cos(phi)*dx + sin(theta)*sin(phi)*dy + cos(theta)*dz; + dtheta = cos(theta)*cos(phi)*dx + cos(theta)*sin(phi)*dy - sin(theta)*dz; + dtheta /= r; + dphi = -sin(phi)*dx + cos(phi)*dy; + dphi /= (r * sin(theta)); + + // Step 3: store conserved quantities + L = r * r * sin(theta) * dphi; + double f = 1.0 - SagA.r_s / r; + double dt_dλ = sqrt((dr*dr)/f + r*r*dtheta*dtheta + r*r*sin(theta)*sin(theta)*dphi*dphi); + E = f * dt_dλ; + } + void step(double dλ, double rs) { + if (r <= rs) return; + rk4Step(*this, dλ, rs); + // convert back to cartesian + this->x = r * sin(theta) * cos(phi); + this->y = r * sin(theta) * sin(phi); + this->z = r * cos(theta); + } +}; + +void raytrace(vector& pixels, int W, int H) { + pixels.resize(W * H * 3); + + // build camera basis + vec3 forward = normalize(camera.target - camera.pos); + vec3 right = normalize(cross(forward, vec3(0,1,0))); + vec3 up = cross(right, forward); + float aspect = float(W) / float(H); + float tanHalfFov = tan(radians(camera.fovY) * 0.5f); + + #pragma omp parallel for schedule(dynamic, 4) + for(int y = 0; y < H; ++y) { + for(int x = 0; x < W; ++x) { + // NDC → screen space in [−1,1] + float u = (2.0f * (x + 0.5f) / float(W) - 1.0f) * aspect * tanHalfFov; + float v = (1.0f - 2.0f * (y + 0.5f) / float(H)) * tanHalfFov; + vec3 dir = normalize(u*right + v*up + forward); + + // construct your Ray + Ray ray(camera.pos, dir); + + const int MAX_STEPS = 10000; + const double D_LAMBDA = 1e7; + const double ESCAPE_R = 1e14; + + // 2) march the ray forward in λ + vec3 color(0.0f); + if (!useGeodesics) { + double b = 2.0 * dot(camera.pos, dir); + double c0 = dot(camera.pos, camera.pos) - SagA.r_s*SagA.r_s; + double disc = b*b - 4.0*c0; + if (disc > 0.0) { + double t1 = (-b - sqrt(disc)) * 0.5; + double t2 = (-b + sqrt(disc)) * 0.5; + if (t1 > 0.0 || t2 > 0.0) + color = vec3(1.0f, 0.0f, 0.0f); + } + } + else { + // full null‐geodesic march + Ray ray(camera.pos, dir); + for(int i = 0; i < MAX_STEPS; ++i) { + if (SagA.Intercept(ray.x, ray.y, ray.z)) { + color = vec3(1.0f, 0.0f, 0.0f); + break; + } + ray.step(D_LAMBDA, SagA.r_s); + if (ray.r > ESCAPE_R) { + // escaped to infinity → remains black + break; + } + } + } + + int idx = (y * W + x) * 3; + pixels[idx+0] = (unsigned char)(color.r * 255); + pixels[idx+1] = (unsigned char)(color.g * 255); + pixels[idx+2] = (unsigned char)(color.b * 255); + } + } +} + +void geodesicRHS(const Ray& ray, double rhs[6], double rs) { + double r = ray.r; + double theta = ray.theta; + double dr = ray.dr; + double dtheta = ray.dtheta; + double dphi = ray.dphi; + double E = ray.E; + + double f = 1.0 - rs / r; + double dt_dlambda = E / f; + + // First derivatives + rhs[0] = dr; + rhs[1] = dtheta; + rhs[2] = dphi; + + // Second derivatives (from 3D Schwarzschild null geodesics): + rhs[3] = + - (rs / (2 * r * r)) * f * dt_dlambda * dt_dlambda + + (rs / (2 * r * r * f)) * dr * dr + + r * (dtheta * dtheta + sin(theta) * sin(theta) * dphi * dphi); + + rhs[4] = + - (2.0 / r) * dr * dtheta + + sin(theta) * cos(theta) * dphi * dphi; + + rhs[5] = + - (2.0 / r) * dr * dphi + - 2.0 * cos(theta) / sin(theta) * dtheta * dphi; +} +void addState(const double a[6], const double b[6], double factor, double out[6]) { + for (int i = 0; i < 6; i++) + out[i] = a[i] + b[i] * factor; +} +void rk4Step(Ray& ray, double dλ, double rs) { + double y0[6] = { ray.r, ray.theta, ray.phi, ray.dr, ray.dtheta, ray.dphi }; + double k1[6], k2[6], k3[6], k4[6], temp[6]; + + geodesicRHS(ray, k1, rs); + addState(y0, k1, dλ/2.0, temp); + Ray r2 = ray; + r2.r = temp[0]; r2.theta = temp[1]; r2.phi = temp[2]; + r2.dr = temp[3]; r2.dtheta = temp[4]; r2.dphi = temp[5]; + geodesicRHS(r2, k2, rs); + + addState(y0, k2, dλ/2.0, temp); + Ray r3 = ray; + r3.r = temp[0]; r3.theta = temp[1]; r3.phi = temp[2]; + r3.dr = temp[3]; r3.dtheta = temp[4]; r3.dphi = temp[5]; + geodesicRHS(r3, k3, rs); + + addState(y0, k3, dλ, temp); + Ray r4 = ray; + r4.r = temp[0]; r4.theta = temp[1]; r4.phi = temp[2]; + r4.dr = temp[3]; r4.dtheta = temp[4]; r4.dphi = temp[5]; + geodesicRHS(r4, k4, rs); + + ray.r += (dλ/6.0)*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); + ray.theta += (dλ/6.0)*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1]); + ray.phi += (dλ/6.0)*(k1[2] + 2*k2[2] + 2*k3[2] + k4[2]); + ray.dr += (dλ/6.0)*(k1[3] + 2*k2[3] + 2*k3[3] + k4[3]); + ray.dtheta += (dλ/6.0)*(k1[4] + 2*k2[4] + 2*k3[4] + k4[4]); + ray.dphi += (dλ/6.0)*(k1[5] + 2*k2[5] + 2*k3[5] + k4[5]); +} + +void setupCameraCallbacks(GLFWwindow* window) { + glfwSetWindowUserPointer(window, &camera); + glfwSetMouseButtonCallback(window, Camera::mouseButtonCallback); + glfwSetCursorPosCallback(window, Camera::cursorPosCallback); + glfwSetScrollCallback(window, Camera::scrollCallback); + glfwSetKeyCallback(window, Engine::keyCallback); +} + +// -- MAIN -- // +int main() { + setupCameraCallbacks(engine.window); + vector pixels(engine.WIDTH * engine.HEIGHT * 3); + + auto t0 = Clock::now(); + lastPrintTime = std::chrono::duration(t0.time_since_epoch()).count(); + + while (!glfwWindowShouldClose(engine.window)) { + raytrace(pixels, engine.WIDTH, engine.HEIGHT); + engine.renderScene(pixels, engine.WIDTH, engine.HEIGHT); + + // 2) FPS counting + framesCount++; + auto t1 = Clock::now(); + double now = std::chrono::duration(t1.time_since_epoch()).count(); + if (now - lastPrintTime >= 1.0) { + cout << "FPS: " << framesCount / (now - lastPrintTime) << "\n"; + framesCount = 0; + lastPrintTime = now; + } + + } + + glfwDestroyWindow(engine.window); + glfwTerminate(); + return 0; +} + + + + + + + + + + + + + + + + + + + // 2) FPS counting + // framesCount++; + // auto t1 = Clock::now(); + // double now = std::chrono::duration(t1.time_since_epoch()).count(); + // if (now - lastPrintTime >= 1.0) { + // cout << "FPS: " << framesCount / (now - lastPrintTime) << "\n"; + // framesCount = 0; + // lastPrintTime = now; + // } + //raytrace(pixels, engine.WIDTH, engine.HEIGHT); diff --git a/Gravity_Sim.zip b/Gravity_Sim.zip new file mode 100644 index 0000000..098c42d Binary files /dev/null and b/Gravity_Sim.zip differ diff --git a/README.md b/README.md index be6ec2c..aae4faf 100644 --- a/README.md +++ b/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 :) diff --git a/black_hole.cpp b/black_hole.cpp index a9708f5..8999365 100644 --- a/black_hole.cpp +++ b/black_hole.cpp @@ -5,57 +5,325 @@ #include #include #include +#define _USE_MATH_DEFINES #include -// #include -// #include -// #include - +#include +#include +#include +#include +#include +#include +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif using namespace glm; using namespace std; +using Clock = std::chrono::high_resolution_clock; -// global vars -const int WIDTH = 800; -const int HEIGHT = 600; -const float G = 6.67430 * pow(10, -11); +// VARS +double lastPrintTime = 0.0; +int framesCount = 0; +double c = 299792458.0; +double G = 6.67430e-11; +struct Ray; +bool Gravity = false; -// functions +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; -// structures and classes :D -class Engine{ -public: - // -- Quad & Texture render + float azimuth = 0.0f; + float elevation = M_PI / 2.0f; + + float orbitSpeed = 0.01f; + float panSpeed = 0.01f; + double zoomSpeed = 25e9f; + + bool dragging = false; + bool panning = false; + bool moving = false; // For compute shader optimization + double lastX = 0.0, lastY = 0.0; + + // Calculate camera position in world space + vec3 position() const { + float clampedElevation = clamp(elevation, 0.01f, float(M_PI) - 0.01f); + // Orbit around (0,0,0) always + return vec3( + radius * sin(clampedElevation) * cos(azimuth), + radius * cos(clampedElevation), + radius * sin(clampedElevation) * sin(azimuth) + ); + } + void update() { + // Always keep target at black hole center + target = vec3(0.0f, 0.0f, 0.0f); + if(dragging | panning) { + moving = true; + } else { + moving = false; + } + } + + void processMouseMove(double x, double y) { + float dx = float(x - lastX); + float dy = float(y - lastY); + + if (dragging && panning) { + // Pan: Shift + Left or Middle Mouse + // Disable panning to keep camera centered on black hole + } + else if (dragging && !panning) { + // Orbit: Left mouse only + azimuth += dx * orbitSpeed; + elevation -= dy * orbitSpeed; + elevation = clamp(elevation, 0.01f, float(M_PI) - 0.01f); + } + + lastX = x; + lastY = y; + update(); + } + void processMouseButton(int button, int action, int mods, GLFWwindow* win) { + if (button == GLFW_MOUSE_BUTTON_LEFT || button == GLFW_MOUSE_BUTTON_MIDDLE) { + if (action == GLFW_PRESS) { + dragging = true; + // Disable panning so camera always orbits center + panning = false; + glfwGetCursorPos(win, &lastX, &lastY); + } else if (action == GLFW_RELEASE) { + dragging = false; + panning = false; + } + } + if (button == GLFW_MOUSE_BUTTON_RIGHT) { + if (action == GLFW_PRESS) { + Gravity = true; + } else if (action == GLFW_RELEASE) { + Gravity = false; + } + } + } + void processScroll(double xoffset, double yoffset) { + radius -= yoffset * zoomSpeed; + radius = clamp(radius, minRadius, maxRadius); + update(); + } + void processKey(int key, int scancode, int action, int mods) { + if (action == GLFW_PRESS && key == GLFW_KEY_G) { + Gravity = !Gravity; + cout << "[INFO] Gravity turned " << (Gravity ? "ON" : "OFF") << endl; + } + } +}; +Camera camera; + +struct BlackHole { + vec3 position; + double mass; + double radius; + double r_s; + + BlackHole(vec3 pos, float m) : position(pos), mass(m) {r_s = 2.0 * G * mass / (c*c);} + bool Intercept(float px, float py, float pz) const { + double dx = double(px) - double(position.x); + double dy = double(py) - double(position.y); + double dz = double(pz) - double(position.z); + double dist2 = dx * dx + dy * dy + dz * dz; + return dist2 < r_s * r_s; + } +}; +BlackHole SagA(vec3(0.0f, 0.0f, 0.0f), 8.54e36); // Sagittarius A black hole +struct ObjectData { + vec4 posRadius; // xyz = position, w = radius + vec4 color; // rgb = color, a = unused + float mass; + vec3 velocity = vec3(0.0f, 0.0f, 0.0f); // Initial velocity +}; +vector objects = { + { vec4(4e11f, 0.0f, 0.0f, 4e10f) , vec4(1,1,0,1), 1.98892e30 }, + { vec4(0.0f, 0.0f, 4e11f, 4e10f) , vec4(1,0,0,1), 1.98892e30 }, + { vec4(0.0f, 0.0f, 0.0f, SagA.r_s) , vec4(0,0,0,1), SagA.mass }, + //{ vec4(6e10f, 0.0f, 0.0f, 5e10f), vec4(0,1,0,1) } +}; + +struct Engine { + GLuint gridShaderProgram; + // -- Quad & Texture render -- // GLFWwindow* window; GLuint quadVAO; GLuint texture; GLuint shaderProgram; + GLuint computeProgram = 0; + // -- UBOs -- // + GLuint cameraUBO = 0; + GLuint diskUBO = 0; + GLuint objectsUBO = 0; + // -- grid mess vars -- // + GLuint gridVAO = 0; + GLuint gridVBO = 0; + GLuint gridEBO = 0; + int gridIndexCount = 0; - Engine(){ - this->window = StartGLFW(); + 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]; } - GLFWwindow* StartGLFW(){ - if(!glfwInit()){ - std::cerr<<"glfw failed init, PANIC PANIC!"<& objects) { + const int gridSize = 25; + const float spacing = 1e10f; // tweak this + + vector vertices; + vector indices; + + for (int z = 0; z <= gridSize; ++z) { + for (int x = 0; x <= gridSize; ++x) { + float worldX = (x - gridSize / 2) * spacing; + float worldZ = (z - gridSize / 2) * spacing; + + float y = 0.0f; + + // ✅ Warp grid using Schwarzschild geometry + for (const auto& obj : objects) { + vec3 objPos = vec3(obj.posRadius); + double mass = obj.mass; + double radius = obj.posRadius.w; + + double r_s = 2.0 * G * mass / (c * c); + double dx = worldX - objPos.x; + double dz = worldZ - objPos.z; + double dist = sqrt(dx * dx + dz * dz); + + // prevent sqrt of negative or divide-by-zero (inside or at the black hole center) + if (dist > r_s) { + double deltaY = 2.0 * sqrt(r_s * (dist - r_s)); + y += static_cast(deltaY) - 3e10f; + } else { + // 🔴 For points inside or at r_s: make it dip down sharply + y += 2.0f * static_cast(sqrt(r_s * r_s)) - 3e10f; // or add a deep pit + } + } + + vertices.emplace_back(worldX, y, worldZ); + } } - 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; + // 🧩 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); + } } - glViewport(0, 0, WIDTH, HEIGHT); - return window; - }; + // 🔌 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 @@ -96,8 +364,198 @@ public: 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(); - std::vector QuadVAO(){ + GLuint shader = glCreateShader(type); + glShaderSource(shader, 1, &src, nullptr); + glCompileShader(shader); + + GLint success; + glGetShaderiv(shader, GL_COMPILE_STATUS, &success); + if (!success) { + GLint logLen; + glGetShaderiv(shader, GL_INFO_LOG_LENGTH, &logLen); + std::vector log(logLen); + glGetShaderInfoLog(shader, logLen, nullptr, log.data()); + std::cerr << "Shader compile error (" << path << "):\n" << log.data() << "\n"; + exit(EXIT_FAILURE); + } + return shader; + }; + + GLuint vertShader = loadShader(vertPath, GL_VERTEX_SHADER); + GLuint fragShader = loadShader(fragPath, GL_FRAGMENT_SHADER); + + GLuint program = glCreateProgram(); + glAttachShader(program, vertShader); + glAttachShader(program, fragShader); + glLinkProgram(program); + + GLint linkSuccess; + glGetProgramiv(program, GL_LINK_STATUS, &linkSuccess); + if (!linkSuccess) { + GLint logLen; + glGetProgramiv(program, GL_INFO_LOG_LENGTH, &logLen); + std::vector log(logLen); + glGetProgramInfoLog(program, logLen, nullptr, log.data()); + std::cerr << "Shader link error:\n" << log.data() << "\n"; + exit(EXIT_FAILURE); + } + + glDeleteShader(vertShader); + glDeleteShader(fragShader); + + return program; + } + GLuint CreateComputeProgram(const char* path) { + // 1) read GLSL source + std::ifstream in(path); + if(!in.is_open()) { + std::cerr << "Failed to open compute shader: " << path << "\n"; + exit(EXIT_FAILURE); + } + std::stringstream ss; + ss << in.rdbuf(); + std::string srcStr = ss.str(); + const char* src = srcStr.c_str(); + + // 2) compile + GLuint cs = glCreateShader(GL_COMPUTE_SHADER); + glShaderSource(cs, 1, &src, nullptr); + glCompileShader(cs); + GLint ok; + glGetShaderiv(cs, GL_COMPILE_STATUS, &ok); + if(!ok) { + GLint logLen; + glGetShaderiv(cs, GL_INFO_LOG_LENGTH, &logLen); + std::vector log(logLen); + glGetShaderInfoLog(cs, logLen, nullptr, log.data()); + std::cerr << "Compute shader compile error:\n" << log.data() << "\n"; + exit(EXIT_FAILURE); + } + + // 3) link + GLuint prog = glCreateProgram(); + glAttachShader(prog, cs); + glLinkProgram(prog); + glGetProgramiv(prog, GL_LINK_STATUS, &ok); + if(!ok) { + GLint logLen; + glGetProgramiv(prog, GL_INFO_LOG_LENGTH, &logLen); + std::vector log(logLen); + glGetProgramInfoLog(prog, logLen, nullptr, log.data()); + std::cerr << "Compute shader link error:\n" << log.data() << "\n"; + exit(EXIT_FAILURE); + } + + glDeleteShader(cs); + return prog; + } + void dispatchCompute(const Camera& cam) { + // determine target compute‐res + int cw = cam.moving ? COMPUTE_WIDTH : 200; + int ch = cam.moving ? COMPUTE_HEIGHT : 150; + + // 1) reallocate the texture if needed + glBindTexture(GL_TEXTURE_2D, texture); + glTexImage2D(GL_TEXTURE_2D, + 0, // mip + GL_RGBA8, // internal format + cw, // width + ch, // height + 0, GL_RGBA, + GL_UNSIGNED_BYTE, + nullptr); + + // 2) bind compute program & UBOs + glUseProgram(computeProgram); + uploadCameraUBO(cam); + uploadDiskUBO(); + uploadObjectsUBO(objects); + + // 3) bind it as image unit 0 + glBindImageTexture(0, texture, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA8); + + // 4) dispatch grid + GLuint groupsX = (GLuint)std::ceil(cw / 16.0f); + GLuint groupsY = (GLuint)std::ceil(ch / 16.0f); + glDispatchCompute(groupsX, groupsY, 1); + + // 5) sync + glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT); + } + void uploadCameraUBO(const Camera& cam) { + struct UBOData { + vec3 pos; float _pad0; + vec3 right; float _pad1; + vec3 up; float _pad2; + vec3 forward; float _pad3; + float tanHalfFov; + float aspect; + bool moving; + int _pad4; + } data; + vec3 fwd = normalize(cam.target - cam.position()); + vec3 up = vec3(0, 1, 0); // y axis is up, so disk is in x-z plane + vec3 right = normalize(cross(fwd, up)); + up = cross(right, fwd); + + data.pos = cam.position(); + data.right = right; + data.up = up; + data.forward = fwd; + data.tanHalfFov = tan(radians(60.0f * 0.5f)); + data.aspect = float(WIDTH) / float(HEIGHT); + data.moving = cam.dragging || cam.panning; + + glBindBuffer(GL_UNIFORM_BUFFER, cameraUBO); + glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(UBOData), &data); + } + void uploadObjectsUBO(const vector& objs) { + struct UBOData { + int numObjects; + float _pad0, _pad1, _pad2; // <-- pad out to 16 bytes + vec4 posRadius[16]; + vec4 color[16]; + float mass[16]; + } data; + + size_t count = std::min(objs.size(), size_t(16)); + data.numObjects = static_cast(count); + + for (size_t i = 0; i < count; ++i) { + data.posRadius[i] = objs[i].posRadius; + data.color[i] = objs[i].color; + data.mass[i] = objs[i].mass; + } + + // Upload + glBindBuffer(GL_UNIFORM_BUFFER, objectsUBO); + glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(data), &data); + } + void uploadDiskUBO() { + // disk + float r1 = SagA.r_s * 2.2f; // inner radius just outside the event horizon + float r2 = SagA.r_s * 5.2f; // outer radius of the disk + float num = 2.0; // number of rays + float thickness = 1e9f; // padding for std140 alignment + float diskData[4] = { r1, r2, num, thickness }; + + glBindBuffer(GL_UNIFORM_BUFFER, diskUBO); + glBufferSubData(GL_UNIFORM_BUFFER, 0, sizeof(diskData), diskData); + } + + vector QuadVAO(){ float quadVertices[] = { // positions // texCoords -1.0f, 1.0f, 0.0f, 1.0f, // top left @@ -107,7 +565,6 @@ public: -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; @@ -128,359 +585,126 @@ public: glBindTexture(GL_TEXTURE_2D, texture); glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); - std::vector VAOtexture = {VAO, texture}; + glBindTexture(GL_TEXTURE_2D, texture); + glTexImage2D(GL_TEXTURE_2D, + 0, // mip + GL_RGBA8, // internal format + COMPUTE_WIDTH, + COMPUTE_HEIGHT, + 0, + GL_RGBA, + GL_UNSIGNED_BYTE, + nullptr); + vector VAOtexture = {VAO, texture}; return VAOtexture; } - void renderScene(const std::vector& pixels, int texWidth, int texHeight) { - // update texture w/ ray-tracing results - glBindTexture(GL_TEXTURE_2D, texture); - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, texWidth, texHeight, 0, GL_RGB, GL_UNSIGNED_BYTE, pixels.data()); - - // clear screen and draw textured quad + void renderScene() { glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); glUseProgram(shaderProgram); - - GLint textureLocation = glGetUniformLocation(shaderProgram, "screenTexture"); - glUniform1i(textureLocation, 0); - 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(); }; - std::vector OptimizeMovement(double lastMovementTime){ - double currentTime = glfwGetTime(); - bool isMoving = (currentTime - lastMovementTime < 0.2); - int renderFactor = isMoving ? 4 : 2; - int rWidth = WIDTH / renderFactor; - int rHeight = HEIGHT / renderFactor; - std::vector vec = {rWidth, rHeight}; - return vec; - } }; -class Camera{ -public: - vec3 target; - float distance; - float pitch; - float yaw; - vec3 position; - vec3 up; - - // mouse handling - bool middleMousePressed = false; - double lastX = 0.0, lastY = 0.0; - float orbitSpeed = 0.4f; - float zoomSpeed = 2.0f; +Engine engine; +void setupCameraCallbacks(GLFWwindow* window) { + glfwSetWindowUserPointer(window, &camera); - float fov = 60.0f; - double lastMovementTime = 0.0; - // default: look at (0,0,0), 5 units far. - Camera(vec3 t = vec3(0.0f, 0.0f, -19.0f), float dist = 5.0f, float yawVal = -90.0f, float pitchVal = 0.0f, float fovVal = 90.0f) - : target(t), distance(dist), yaw(yawVal), pitch(pitchVal), fov(fovVal) { - up = vec3(0, 1, 0); - updatePosition(); - } - void updatePosition() { - float radYaw = radians(yaw); - float radPitch = radians(pitch); - position.x = target.x + distance * cos(radPitch) * cos(radYaw); - position.y = target.y + distance * sin(radPitch); - position.z = target.z + distance * cos(radPitch) * sin(radYaw); - } - - // Member function to handle mouse button events. - void handleMouseButton(int button, int action, int mods, GLFWwindow* window) { - if (button == GLFW_MOUSE_BUTTON_MIDDLE) { - if (action == GLFW_PRESS) { - middleMousePressed = true; - glfwGetCursorPos(window, &lastX, &lastY); - lastMovementTime = glfwGetTime(); - } else if (action == GLFW_RELEASE) { - middleMousePressed = false; - } - } - } - void handleCursorPosition(double xpos, double ypos, GLFWwindow* window) { - if (!middleMousePressed) - return; + glfwSetMouseButtonCallback(window, [](GLFWwindow* win, int button, int action, int mods) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(win); + cam->processMouseButton(button, action, mods, win); + }); - double deltaX = xpos - lastX; - double deltaY = -(ypos - lastY); + glfwSetCursorPosCallback(window, [](GLFWwindow* win, double x, double y) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(win); + cam->processMouseMove(x, y); + }); - // If shift is held, pan the camera's target. - if (glfwGetKey(window, GLFW_KEY_LEFT_SHIFT) == GLFW_PRESS || - glfwGetKey(window, GLFW_KEY_RIGHT_SHIFT) == GLFW_PRESS) { - vec3 forward = normalize(target - position); - vec3 right = normalize(cross(forward, up)); - vec3 camUp = cross(right, forward); - float panSpeed = 0.005f * distance; - target += -right * (float)deltaX * panSpeed + camUp * (float)deltaY * panSpeed; - } - // Otherwise, orbit the camera. - else { - yaw += (float)deltaX * orbitSpeed; - pitch += (float)deltaY * orbitSpeed; - if (pitch > 89.0f) pitch = 89.0f; - if (pitch < -89.0f) pitch = -89.0f; - } - updatePosition(); - lastX = xpos; - lastY = ypos; - lastMovementTime = glfwGetTime(); - } - void handleScroll(double xoffset, double yoffset, GLFWwindow* window) { - // If this is the first input, initialize mouse position - if (lastX == 0 && lastY == 0) { - glfwGetCursorPos(window, &lastX, &lastY); - } + glfwSetScrollCallback(window, [](GLFWwindow* win, double xoffset, double yoffset) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(win); + cam->processScroll(xoffset, yoffset); + }); - distance -= (float)yoffset * zoomSpeed; - if (distance < 1.0f) - distance = 1.0f; - - updatePosition(); - lastMovementTime = glfwGetTime(); - } - static void mouseButtonCallback(GLFWwindow* window, int button, int action, int mods) { - Camera* cam = static_cast(glfwGetWindowUserPointer(window)); - cam->handleMouseButton(button, action, mods, window); - } - static void cursorPositionCallback(GLFWwindow* window, double xpos, double ypos) { - Camera* cam = static_cast(glfwGetWindowUserPointer(window)); - cam->handleCursorPosition(xpos, ypos, window); - } - static void scrollCallback(GLFWwindow* window, double xoffset, double yoffset) { - Camera* cam = static_cast(glfwGetWindowUserPointer(window)); - cam->handleScroll(xoffset, yoffset, window); - } - - void registerCallbacks(GLFWwindow* window) { - glfwSetWindowUserPointer(window, this); - glfwSetMouseButtonCallback(window, Camera::mouseButtonCallback); - glfwSetCursorPosCallback(window, Camera::cursorPositionCallback); - glfwSetScrollCallback(window, Camera::scrollCallback); - } -}; - -struct Ray{ - vec3 direction; - vec3 origin; - Ray(vec3 o, vec3 d) : origin(o), direction(normalize(d)){} -}; -struct Material{ - vec3 color; - float specular; - float emission; - float shiny; - Material(vec3 c, float s, float e, float sh = 4.0) : color(c), specular(s), emission(e), shiny(sh) {} -}; -struct Object{ - vec3 position; - vec3 velocity; - float radius; - float mass = 7.3 * pow(10, 22); - Material material; - - Object(vec3 p, float r, Material m, vec3 v = vec3(0.0)) : position(p), radius(r), material(m), velocity(v) {} - // raytracing - bool Intersect(Ray &ray, float &t) const { - vec3 oc = ray.origin - position; // centre to origin - float a = dot(ray.direction, ray.direction); // ray straightness / magnitute - float b = 2.0f * dot(oc, ray.direction); // orientation towards - float c = 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 - position); - } - void UpdatePos(){ - this->position[0] += this->velocity[0]; - this->position[1] += this->velocity[1]; - this->position[2] += this->velocity[2]; - } - // gravity - void accelerate(float x, float y, float z){ - this->velocity[0] += x; - this->velocity[1] += y; - this->velocity[2] += z; - } -}; - -class Scene { -public: - std::vector objs; - //std::vector lightPos; - vector lights; - - vec3 trace(Ray &ray, int depth = 0){ - const int maxDepth = 4; - if (depth >= maxDepth) return vec3(0.0f, 0.0f, 0.0f); - - 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){return vec3(0.05f, 0.05f, 0.1f);} - - // hit - vec3 hitPoint = ray.origin + ray.direction * closest; // point on obj hit by ray - vec3 normal = hitObj->getNormal(hitPoint); - vec3 viewDir = normalize(-ray.direction); // direction to camera - vec3 finalColor = hitObj->material.color * 0.3f; - - if (hitObj->material.emission > 0.0f) { - finalColor += hitObj->material.color * hitObj->material.emission; - } - - for(auto& light : lights) { - if (light == hitObj) continue; - vec3 lightDir = normalize(light->position - hitPoint); // light to hitpoint dir - float distanceToLight = length(light->position - hitPoint); - if (distanceToLight < 0.001f) continue; - - // Compute diffuse lighting - vec3 lightColor = light->material.color * light->material.emission; - float diff = std::max(dot(normal, lightDir), 0.0f); - float attenuation = 1.0f / (distanceToLight * distanceToLight); - vec3 diffuse = hitObj->material.color * lightColor * diff * attenuation; - - // Compute specular Lighting - vec3 reflectDir = reflect(-lightDir, normal); - float spec = pow(std::max(dot(viewDir, reflectDir), 0.0f), hitObj->material.shiny); - vec3 specular = lightColor * hitObj->material.specular * spec * attenuation; - - Ray shadowRay(hitPoint + normal * 0.001f, lightDir); - bool inShadow = false; - for(const auto& obj : objs) { - if (&obj != hitObj) { - float t; - if (obj.Intersect(shadowRay, t) && t < distanceToLight) { - inShadow = true; - break; - } - } - } - // Add light contribution if not in shadow - if (!inShadow) { - finalColor += diffuse + specular; - } - - } - if(hitObj->material.specular > 0.0f) { - vec3 reflectDir = reflect(ray.direction, normal); // Reflect ray direction - Ray reflectRay(hitPoint + normal * 0.001f, reflectDir); - vec3 reflectedColor = trace(reflectRay, depth + 1); // Recurse - finalColor += reflectedColor * hitObj->material.specular * 0.3f; // Add reflection scaled by specular intensity - } - return finalColor; - }; -}; - -// --- main loop ---- // -int main(){ - Engine engine; - Scene scene; - Camera camera(vec3(0.0f, 0.0f, -9.0f), -15.0f, -90.0f, 0.0f, 90.0f); - camera.registerCallbacks(engine.window); - - scene.objs = { - // position radius material: color specular emission - Object(vec3(0.0f, -5.0f, -19.0f), 12.0f, Material(vec3(1.0f, 0.0f, 0.0f), 0.9f, 10.0f)), - Object(vec3(20.0f, -2.0f, -11.0f), 5.5f, Material(vec3(0.1f, 1.0f, 0.1f), 0.2f, 0.5f)), - Object(vec3(-17.0f, -1.0f, -6.0f), 7.0f, Material(vec3(0.1f, 0.1f, 1.0f), 0.8f, 0.3f)), - Object(vec3(-17.0f, -10.0f, 9.0f), 7.0f, Material(vec3(0.1f, 1.0f, 1.0f), 0.0f, 0.3f)), - }; - // -- loop -- // - double lastFrame = glfwGetTime(); - while(!glfwWindowShouldClose(engine.window)){ - glClear(GL_COLOR_BUFFER_BIT); - double currentTime = glfwGetTime(); - double deltaTime = currentTime - lastFrame; - lastFrame = currentTime; - - int rWidth = engine.OptimizeMovement(camera.lastMovementTime)[0]; - int rHeight = engine.OptimizeMovement(camera.lastMovementTime)[1]; - std::vector pixels(rWidth * rHeight * 3); - - // Update light sources - scene.lights.clear(); - for (const auto& obj : scene.objs) { - if (obj.material.emission > 0.0f) { - scene.lights.push_back(&obj); - } - } - - // render texture (pxl by pxl) - for(int y = 0; y < rHeight; ++y){ - for(int x = 0; x < rWidth; ++x){ - float scale = tan(radians(camera.fov * 0.5f)); - - float aspectRatio = float(rWidth) / float(rHeight); - float u = float(x) / float(rWidth); // % of width - float v = float(y) / float(rHeight); // % of height - - // Convert screen coordinates to camera space coordinates with FOV adjustment - float x_camera = (2.0f * u - 1.0f) * aspectRatio * scale; - float y_camera = (1.0f - 2.0f * v) * scale; // (1 - 2*v) is equivalent to -(2*v - 1) - - // Transform the ray from camera space to world space. - vec3 forward = normalize(camera.target - camera.position); - vec3 right = normalize(cross(forward, vec3(0.0f, 1.0f, 0.0f))); - vec3 up = cross(right, forward); - - vec3 direction = normalize(x_camera * right + y_camera * up + forward); - - Ray ray(camera.position, direction); - vec3 color = scene.trace(ray); - color = color / (color + vec3(0.5f)); // Reinhard tone mapping - color = clamp(color, 0.0f, 1.0f); - - int index = (y * rWidth + x) * 3; - pixels[index + 0] = static_cast(color.r * 255); - pixels[index + 1] = static_cast(color.g * 255); - pixels[index + 2] = static_cast(color.b * 255); - } - } - - for(auto& obj : scene.objs) { - if(obj.position[1] > 0){ - obj.velocity *= -0.8f; - obj.position[1] = 0.1f; - } - obj.accelerate(0.0, 9.81 * deltaTime, 0.0); - obj.UpdatePos(); - } - - engine.renderScene(pixels, rWidth, rHeight); - } - - glfwTerminate(); + glfwSetKeyCallback(window, [](GLFWwindow* win, int key, int scancode, int action, int mods) { + Camera* cam = (Camera*)glfwGetWindowUserPointer(win); + cam->processKey(key, scancode, action, mods); + }); } -// func dec's +// -- MAIN -- // +int main() { + setupCameraCallbacks(engine.window); + vector pixels(engine.WIDTH * engine.HEIGHT * 3); + + auto t0 = Clock::now(); + lastPrintTime = chrono::duration(t0.time_since_epoch()).count(); + + double lastTime = glfwGetTime(); + int renderW = 800, renderH = 600, numSteps = 80000; + while (!glfwWindowShouldClose(engine.window)) { + glClearColor(0.0f, 0.0f, 0.0f, 1.0f); // optional, but good practice + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + + double now = glfwGetTime(); + double dt = now - lastTime; // seconds since last frame + lastTime = now; + + // Gravity + for (auto& obj : objects) { + for (auto& obj2 : objects) { + if (&obj == &obj2) continue; // skip self-interaction + float dx = obj2.posRadius.x - obj.posRadius.x; + float dy = obj2.posRadius.y - obj.posRadius.y; + float dz = obj2.posRadius.z - obj.posRadius.z; + float distance = sqrt(dx * dx + dy * dy + dz * dz); + if (distance > 0) { + vector direction = {dx / distance, dy / distance, dz / distance}; + //distance *= 1000; + double Gforce = (G * obj.mass * obj2.mass) / (distance * distance); + + double acc1 = Gforce / obj.mass; + std::vector acc = {direction[0] * acc1, direction[1] * acc1, direction[2] * acc1}; + if (Gravity) { + obj.velocity.x += acc[0]; + obj.velocity.y += acc[1]; + obj.velocity.z += acc[2]; + + obj.posRadius.x += obj.velocity.x; + obj.posRadius.y += obj.velocity.y; + obj.posRadius.z += obj.velocity.z; + cout << "velocity: " <= disk_r1 && r <= disk_r2); +} + +void main() { + int WIDTH = cam.moving ? 200 : 200; + int HEIGHT = cam.moving ? 150 : 150; + + ivec2 pix = ivec2(gl_GlobalInvocationID.xy); + if (pix.x >= WIDTH || pix.y >= HEIGHT) return; + + // Init Ray + float u = (2.0 * (pix.x + 0.5) / WIDTH - 1.0) * cam.aspect * cam.tanHalfFov; + float v = (1.0 - 2.0 * (pix.y + 0.5) / HEIGHT) * cam.tanHalfFov; + vec3 dir = normalize(u * cam.camRight - v * cam.camUp + cam.camForward); + Ray ray = initRay(cam.camPos, dir); + + vec4 color = vec4(0.0); + vec3 prevPos = vec3(ray.x, ray.y, ray.z); + float lambda = 0.0; + + bool hitBlackHole = false; + bool hitDisk = false; + bool hitObject = false; + + int steps = cam.moving ? 60000 : 60000; + + for (int i = 0; i < steps; ++i) { + if (intercept(ray, SagA_rs)) { hitBlackHole = true; break; } + rk4Step(ray, D_LAMBDA); + lambda += D_LAMBDA; + + vec3 newPos = vec3(ray.x, ray.y, ray.z); + if (crossesEquatorialPlane(prevPos, newPos)) { hitDisk = true; break; } + if (interceptObject(ray)) { hitObject = true; break; } + prevPos = newPos; + if (ray.r > ESCAPE_R) break; + } + + if (hitDisk) { + double r = length(vec3(ray.x, ray.y, ray.z)) / disk_r2; + vec3 diskColor = vec3(1.0, r, 0.2); + //r = 1.0 - abs(r - 0.5) * 2.0; + color = vec4(diskColor, r); + + } else if (hitBlackHole) { + color = vec4(0.0, 0.0, 0.0, 1.0); + + } else if (hitObject) { + // Compute shading + vec3 P = vec3(ray.x, ray.y, ray.z); + vec3 N = normalize(P - hitCenter); + vec3 V = normalize(cam.camPos - P); + float ambient = 0.1; + float diff = max(dot(N, V), 0.0); + float intensity = ambient + (1.0 - ambient) * diff; + vec3 shaded = objectColor.rgb * intensity; + color = vec4(shaded, objectColor.a); + + } else { + color = vec4(0.0); + } + + imageStore(outImage, pix, color); +} diff --git a/vs_code/c_cpp_properties.json b/vs_code/c_cpp_properties.json new file mode 100644 index 0000000..13ff85a --- /dev/null +++ b/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 +} \ No newline at end of file diff --git a/vs_code/launch.json b/vs_code/launch.json new file mode 100644 index 0000000..4d1a665 --- /dev/null +++ b/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" + } + ] +} \ No newline at end of file diff --git a/vs_code/settings.json b/vs_code/settings.json new file mode 100644 index 0000000..145c23d --- /dev/null +++ b/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" + } +} \ No newline at end of file diff --git a/vs_code/tasks.json b/vs_code/tasks.json new file mode 100644 index 0000000..b9a6a36 --- /dev/null +++ b/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" + } + ] +} \ No newline at end of file