diff --git a/.idea/vcs.xml b/.idea/vcs.xml index 94a25f7..35eb1dd 100644 --- a/.idea/vcs.xml +++ b/.idea/vcs.xml @@ -1,6 +1,6 @@ - + \ No newline at end of file diff --git a/CelestialBody.cpp b/CelestialBody.cpp index 1919b9e..c6ab7f0 100644 --- a/CelestialBody.cpp +++ b/CelestialBody.cpp @@ -12,10 +12,10 @@ std::string vec3_to_string(const glm::vec3& v) { return ss.str(); } -CelestialBody::CelestialBody(float mass, const glm::vec3& position, const glm::vec3& velocity) +CelestialBody::CelestialBody(double mass, const glm::dvec3& position, const glm::dvec3& velocity) : mass(mass), position(position), velocity(velocity), acceleration(0.0f) {} -void CelestialBody::update(float dt) { +void CelestialBody::update(double dt) { if (glm::any(glm::isnan(velocity)) || glm::any(glm::isinf(velocity))) { std::cout << "Warning: Invalid velocity detected: " << vec3_to_string(velocity) << std::endl; velocity = glm::vec3(0.0f); @@ -32,6 +32,13 @@ void CelestialBody::update(float dt) { acceleration = glm::vec3(0.0f); } -void CelestialBody::applyForce(const glm::vec3& force) { +void CelestialBody::applyForce(const glm::dvec3& force) { acceleration += force / mass; +} + +void CelestialBody::addToTrajectory(const glm::dvec3& position) { + trajectory.push_back(position); + if (trajectory.size() > MAX_TRAJECTORY_POINTS) { + trajectory.erase(trajectory.begin()); + } } \ No newline at end of file diff --git a/CelestialBody.h b/CelestialBody.h index 61c31b0..1c22a60 100644 --- a/CelestialBody.h +++ b/CelestialBody.h @@ -7,22 +7,27 @@ #pragma once #include #include +#include class CelestialBody { public: - CelestialBody(float mass, const glm::vec3& position, const glm::vec3& velocity); + CelestialBody(double mass, const glm::dvec3& position, const glm::dvec3& velocity); - void update(float dt); - void applyForce(const glm::vec3& force); + void update(double dt); + void applyForce(const glm::dvec3& force); - float getMass() const { return mass; } - glm::vec3 getPosition() const { return position; } - glm::vec3 getVelocity() const { return velocity; } + [[nodiscard]] double getMass() const { return mass; } + [[nodiscard]] glm::dvec3 getPosition() const { return position; } + [[nodiscard]] glm::dvec3 getVelocity() const { return velocity; } + void addToTrajectory(const glm::dvec3& position); + const std::vector& getTrajectory() const { return trajectory; } private: - float mass; - glm::vec3 position; - glm::vec3 velocity; - glm::vec3 acceleration; + double mass; + glm::dvec3 position; + glm::dvec3 velocity; + glm::dvec3 acceleration; + std::vector trajectory; + static const size_t MAX_TRAJECTORY_POINTS = 1000; }; #endif //GRAVITY_CELESTIALBODY_H diff --git a/README.md b/README.md new file mode 100644 index 0000000..1bbba71 --- /dev/null +++ b/README.md @@ -0,0 +1,55 @@ +# 3D Gravity Simulator Documentation + +## Overview + +This 3D Gravity Simulator is a C++ program that visualizes the gravitational interactions between celestial bodies in a simplified solar system model. It uses OpenGL for rendering and GLFW for window management and user input. + +## Program Structure + +The simulator consists of several key components: + +1. `Simulator`: Handles the physics calculations and updates the positions of celestial bodies. +2. `Renderer`: Manages the 3D rendering of the celestial bodies, trajectories, and grid. +3. `CelestialBody`: Represents individual celestial bodies with properties like mass, position, and velocity. + +## Physics Implementation + +### Gravitational Force + +The simulator uses Newton's law of universal gravitation to calculate the forces between celestial bodies. The gravitational force between two bodies is given by: + +$$ F = G \frac{m_1 m_2}{r^2} $$ + +Where: +- $F$ is the gravitational force between the two bodies +- $G$ is the gravitational constant ($$6.67430 \times 10^{-11} \, \text{N} \cdot \text{m}^2 / \text{kg}^2$$) +- $m_1$ and $m_2$ are the masses of the two bodies +- $r$ is the distance between the centers of the two bodies + +### Motion Update + +The motion of each celestial body is updated using numerical integration. We use a simple Euler method for updating positions and velocities: + +1. Calculate the net force on each body +2. Calculate acceleration: $$ \vec{a} = \frac{\vec{F}}{m} $$ +3. Update velocity: $$ \vec{v}_{new} = \vec{v}_{old} + \vec{a} \Delta t $$ +4. Update position: $$ \vec{x}_{new} = \vec{x}_{old} + \vec{v}_{new} \Delta t $$ + +Where $\Delta t$ is the time step of the simulation. + +## Rendering + +The program uses OpenGL to render the 3D scene: + +- Celestial bodies are represented as spheres with sizes proportional to their masses (using a logarithmic scale). +- A grid is drawn to provide a reference plane. +- Trajectories of the bodies are drawn as lines, fading out over time. +- The camera can be controlled using WASD keys for movement and the mouse for orientation. + +## Limitations and Simplifications + +1. The simulation uses a fixed time step, which can lead to inaccuracies in long-term simulations. +2. The Euler method for numerical integration is simple but can accumulate errors over time. +3. The scale of the celestial bodies and their distances are not to true scale to make visualization easier. +4. Relativistic effects are not considered; the simulation uses classical Newtonian mechanics. + diff --git a/Renderer.cpp b/Renderer.cpp index e429771..e4a15b2 100644 --- a/Renderer.cpp +++ b/Renderer.cpp @@ -9,7 +9,18 @@ #include #include -Renderer::Renderer(int width, int height) { +Renderer::Renderer(int width, int height) + : cameraPos(3e11f, 2e11f, 3e11f), + cameraFront(glm::normalize(glm::vec3(0.0f) - glm::vec3(3e11f, 2e11f, 3e11f))), + cameraUp(0.0f, 1.0f, 0.0f), + cameraSpeed(1e9f), // Reduced speed + mouseSensitivity(0.05f), // Reduced sensitivity + yaw(-45.0f), + pitch(-30.0f), + firstMouse(true), + lastX(width / 2.0f), + lastY(height / 2.0f) +{ if (!glfwInit()) { throw std::runtime_error("Failed to initialize GLFW"); } @@ -32,6 +43,13 @@ Renderer::Renderer(int width, int height) { glEnable(GL_COLOR_MATERIAL); createSphereMesh(1.0f, 20, 20); + + // Set up camera + glfwSetInputMode(window, GLFW_CURSOR, GLFW_CURSOR_DISABLED); + glfwSetWindowUserPointer(window, this); + glfwSetCursorPosCallback(window, [](GLFWwindow* window, double xpos, double ypos) { + static_cast(glfwGetWindowUserPointer(window))->cursorPosCallback(xpos, ypos); + }); } Renderer::~Renderer() { @@ -48,21 +66,58 @@ void Renderer::render(const Simulator& simulator) { glMatrixMode(GL_PROJECTION); glLoadIdentity(); - gluPerspective(45.0, 1024.0 / 768.0, 1e8, 1e12); + gluPerspective(45.0, 1600.0 / 1200.0, 1e9, 1e13); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); - gluLookAt(3e11, 2e11, 3e11, 0, 0, 0, 0, 1, 0); + glm::vec3 center = glm::vec3(0, 0, 0); // Look at the center of the system + gluLookAt(cameraPos.x, cameraPos.y, cameraPos.z, + center.x, center.y, center.z, + cameraUp.x, cameraUp.y, cameraUp.z); + drawGrid(simulator); - drawGrid(); + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + drawTrajectories(simulator.getBodies()); + glDisable(GL_BLEND); const auto& bodies = simulator.getBodies(); + double maxMass = 0; + double minMass = std::numeric_limits::max(); + + // Find the maximum and minimum masses + for (const auto& body : bodies) { + maxMass = std::max(maxMass, body.getMass()); + minMass = std::min(minMass, body.getMass()); + } + + std::cout << "Camera position: " << cameraPos.x << ", " << cameraPos.y << ", " << cameraPos.z << std::endl; + std::cout << "Camera front: " << cameraFront.x << ", " << cameraFront.y << ", " << cameraFront.z << std::endl; + for (size_t i = 0; i < bodies.size(); ++i) { const auto& body = bodies[i]; - float minSize = 2e9f; - float scaleFactor = std::max(std::cbrt(body.getMass()) * 1e-9f, minSize); + glm::dvec3 pos = body.getPosition(); + std::cout << "Body " << i << " position: " << pos.x << ", " << pos.y << ", " << pos.z << std::endl; + } + + // Calculate the log range + double logMinMass = std::log10(minMass); + double logMaxMass = std::log10(maxMass); + double logRange = logMaxMass - logMinMass; + + for (size_t i = 0; i < bodies.size(); ++i) { + const auto& body = bodies[i]; + + // Calculate the scale factor based on mass + double logMass = std::log10(body.getMass()); + double normalizedLogMass = (logMass - logMinMass) / logRange; + float minScale = 5e9f; // Minimum scale to ensure visibility + float maxScale = 5e10f; // Maximum scale to prevent overly large objects + float scaleFactor = minScale + static_cast(normalizedLogMass) * (maxScale - minScale); + + glm::dvec3 pos = body.getPosition(); + glm::vec3 renderPos(static_cast(pos.x), static_cast(pos.y), static_cast(pos.z)); - glm::vec3 pos = body.getPosition(); std::cout << "Rendering body " << i << " ("; switch(i) { case 0: std::cout << "Sun"; break; @@ -86,7 +141,7 @@ void Renderer::render(const Simulator& simulator) { default: glColor3f(1.0f, 1.0f, 1.0f); break; // White for any additional bodies } - drawSphere(body.getPosition(), scaleFactor); + drawSphere(renderPos, scaleFactor); } } @@ -191,6 +246,7 @@ void Renderer::drawDebugTriangle() { glMatrixMode(GL_MODELVIEW); glLoadIdentity(); + gluLookAt(4e11, 3e11, 4e11, 0, 0, 0, 0, 1, 0); glBegin(GL_TRIANGLES); glColor3f(1.0f, 0.0f, 0.0f); @@ -202,14 +258,100 @@ void Renderer::drawDebugTriangle() { glEnd(); } -void Renderer::drawGrid() { +float Renderer::calculateGravityFieldStrength(const glm::vec3& point, const std::vector& bodies) { + float fieldStrength = 0.0f; + const float G = 6.67430e-11f; // Gravitational constant + const float scalingFactor = 1e20f; // Greatly increased scaling factor + for (const auto& body : bodies) { + glm::dvec3 bodyPos = body.getPosition(); + float distance = glm::length(glm::vec3(bodyPos) - point); + if (distance < 1e9f) distance = 1e9f; // Prevent division by zero + fieldStrength += scalingFactor * G * static_cast(body.getMass()) / (distance * distance); + } + return fieldStrength; +} + +void Renderer::drawGrid(const Simulator& simulator) { + const float gridSize = 5e11f; + const int gridLines = 20; + const float lineSpacing = gridSize / gridLines; + glBegin(GL_LINES); - glColor3f(0.2f, 0.2f, 0.2f); // Gray color for the grid - for (float i = -5e11f; i <= 5e11f; i += 5e10f) { - glVertex3f(i, 0, -5e11f); - glVertex3f(i, 0, 5e11f); - glVertex3f(-5e11f, 0, i); - glVertex3f(5e11f, 0, i); + glColor3f(0.2f, 0.2f, 0.2f); // Lighter gray for better visibility + + for (int i = -gridLines/2; i <= gridLines/2; ++i) { + float pos = i * lineSpacing; + glVertex3f(-gridSize/2, 0, pos); + glVertex3f(gridSize/2, 0, pos); + glVertex3f(pos, 0, -gridSize/2); + glVertex3f(pos, 0, gridSize/2); + } + + glEnd(); +} + +void Renderer::drawTrajectories(const std::vector& bodies) { + glBegin(GL_LINES); + for (const auto& body : bodies) { + const auto& trajectory = body.getTrajectory(); + if (trajectory.size() < 2) continue; + + for (size_t i = 1; i < trajectory.size(); ++i) { + glm::vec3 p1(trajectory[i-1]); + glm::vec3 p2(trajectory[i]); + + // Fade out older parts of the trajectory + float alpha = static_cast(i) / trajectory.size(); + glColor4f(1.0f, 1.0f, 1.0f, alpha * 0.5f); + + glVertex3f(p1.x, p1.y, p1.z); + glVertex3f(p2.x, p2.y, p2.z); + } } glEnd(); +} + +void Renderer::processInput() { + if (glfwGetKey(window, GLFW_KEY_W) == GLFW_PRESS) + cameraPos += cameraSpeed * cameraFront; + if (glfwGetKey(window, GLFW_KEY_S) == GLFW_PRESS) + cameraPos -= cameraSpeed * cameraFront; + if (glfwGetKey(window, GLFW_KEY_A) == GLFW_PRESS) + cameraPos -= glm::normalize(glm::cross(cameraFront, cameraUp)) * cameraSpeed; + if (glfwGetKey(window, GLFW_KEY_D) == GLFW_PRESS) + cameraPos += glm::normalize(glm::cross(cameraFront, cameraUp)) * cameraSpeed; +} + +void Renderer::cursorPosCallback(double xpos, double ypos) { + if (firstMouse) { + lastX = xpos; + lastY = ypos; + firstMouse = false; + } + + float xoffset = xpos - lastX; + float yoffset = lastY - ypos; + lastX = xpos; + lastY = ypos; + + xoffset *= mouseSensitivity; + yoffset *= mouseSensitivity; + + yaw += xoffset; + pitch += yoffset; + + if (pitch > 89.0f) + pitch = 89.0f; + if (pitch < -89.0f) + pitch = -89.0f; + + updateCameraVectors(); +} + +void Renderer::updateCameraVectors() { + glm::vec3 front; + front.x = cos(glm::radians(yaw)) * cos(glm::radians(pitch)); + front.y = sin(glm::radians(pitch)); + front.z = sin(glm::radians(yaw)) * cos(glm::radians(pitch)); + cameraFront = glm::normalize(front); } \ No newline at end of file diff --git a/Renderer.h b/Renderer.h index 2f18a10..c0b02d1 100644 --- a/Renderer.h +++ b/Renderer.h @@ -18,6 +18,8 @@ public: void render(const Simulator& simulator); bool shouldClose(); void swapBuffers(); + void processInput(); + void cursorPosCallback(double xpos, double ypos); private: GLFWwindow* window; @@ -28,6 +30,21 @@ private: GLuint sphereVAO, sphereVBO, sphereEBO; int sphereVertexCount, sphereIndexCount; - void drawGrid(); + void drawGrid(const Simulator& simulator); + float calculateGravityFieldStrength(const glm::vec3& point, const std::vector& bodies); + void drawGravityField(const Simulator& simulator); + void drawTrajectories(const std::vector& bodies); + + glm::vec3 cameraPos; + glm::vec3 cameraFront; + glm::vec3 cameraUp; + float cameraSpeed; + float mouseSensitivity; + float yaw; + float pitch; + bool firstMouse; + double lastX, lastY; + + void updateCameraVectors(); }; #endif //GRAVITY_RENDERER_H diff --git a/Simulator.cpp b/Simulator.cpp index db6204a..f439ff7 100644 --- a/Simulator.cpp +++ b/Simulator.cpp @@ -4,6 +4,7 @@ #include "Simulator.h" #include #include +#include Simulator::Simulator() {} @@ -11,7 +12,12 @@ void Simulator::addBody(const CelestialBody& body) { bodies.push_back(body); } -void Simulator::update(float dt) { +void Simulator::update(double dt) { + // Sort bodies by mass (descending order) + std::sort(bodies.begin(), bodies.end(), [](const CelestialBody& a, const CelestialBody& b) { + return a.getMass() > b.getMass(); + }); + // Calculate and apply gravitational forces for (size_t i = 0; i < bodies.size(); ++i) { glm::vec3 totalForce(0.0f); @@ -25,14 +31,15 @@ void Simulator::update(float dt) { } // Update positions and velocities - for (auto& body : bodies) { - body.update(dt); + for (size_t i = 1; i < bodies.size(); ++i) { // Start from 1 to skip the Sun + bodies[i].update(dt); + bodies[i].addToTrajectory(bodies[i].getPosition()); } } -glm::vec3 Simulator::calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2) { - glm::vec3 direction = body2.getPosition() - body1.getPosition(); - float distance = glm::length(direction); +glm::dvec3 Simulator::calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2) { + glm::dvec3 direction = body2.getPosition() - body1.getPosition(); + double distance = glm::length(direction); // Avoid division by zero and unrealistic forces at very small distances if (distance < 1e9) { @@ -41,13 +48,13 @@ glm::vec3 Simulator::calculateGravitationalForce(const CelestialBody& body1, con } // Use the actual G value - const float G = 6.67430e-11f; - float forceMagnitude = G * (body1.getMass() * body2.getMass()) / (distance * distance); + const double G = 6.67430e-11; + double forceMagnitude = G * (body1.getMass() * body2.getMass()) / (distance * distance); if (std::isnan(forceMagnitude) || std::isinf(forceMagnitude)) { std::cout << "Warning: Invalid force magnitude calculated. Distance: " << distance << ", Masses: " << body1.getMass() << ", " << body2.getMass() << std::endl; - return glm::vec3(0.0f); + return glm::dvec3(0.0); } return glm::normalize(direction) * forceMagnitude; diff --git a/Simulator.h b/Simulator.h index 016b73f..8d40590 100644 --- a/Simulator.h +++ b/Simulator.h @@ -13,13 +13,13 @@ public: Simulator(); void addBody(const CelestialBody& body); - void update(float dt); + void update(double dt); const std::vector& getBodies() const { return bodies; } + glm::dvec3 calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2); + private: std::vector bodies; const float G = 6.67430e-11f; // Gravitational constant - - glm::vec3 calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2); }; #endif //GRAVITY_SIMULATOR_H diff --git a/image.png b/image.png new file mode 100644 index 0000000..80a92b9 Binary files /dev/null and b/image.png differ diff --git a/main.cpp b/main.cpp index bc04310..fde4020 100644 --- a/main.cpp +++ b/main.cpp @@ -8,31 +8,32 @@ int main() { Simulator simulator; - Renderer renderer(1024, 768); // Increased window size for better visibility + Renderer renderer(1600, 1200); // Increased window size for better visibility - // Sun - simulator.addBody(CelestialBody(1.989e30f, glm::vec3(0, 0, 0), glm::vec3(0, 0, 0))); +// Sun (at the center) + simulator.addBody(CelestialBody(1.989e30f, glm::dvec3(0, 0, 0), glm::dvec3(0, 0, 0))); // Mercury - simulator.addBody(CelestialBody(3.285e23f, glm::vec3(57.9e9f, 0, 0), glm::vec3(0, 47.36e3f, 0))); + simulator.addBody(CelestialBody(3.285e23f, glm::dvec3(57.9e9f, 0, 0), glm::dvec3(0, 47.36e3f, 0))); // Venus - simulator.addBody(CelestialBody(4.867e24f, glm::vec3(108.2e9f, 0, 0), glm::vec3(0, 35.02e3f, 0))); + simulator.addBody(CelestialBody(4.867e24f, glm::dvec3(108.2e9f, 0, 0), glm::dvec3(0, 35.02e3f, 0))); // Earth - simulator.addBody(CelestialBody(5.972e24f, glm::vec3(149.6e9f, 0, 0), glm::vec3(0, 29.78e3f, 0))); + simulator.addBody(CelestialBody(5.972e24f, glm::dvec3(149.6e9f, 0, 0), glm::dvec3(0, 29.78e3f, 0))); // Mars - simulator.addBody(CelestialBody(6.39e23f, glm::vec3(227.9e9f, 0, 0), glm::vec3(0, 24.07e3f, 0))); + simulator.addBody(CelestialBody(6.39e23f, glm::dvec3(227.9e9f, 0, 0), glm::dvec3(0, 24.07e3f, 0))); const float dt = 3600.0f; // Time step of 1 hour while (!renderer.shouldClose()) { + renderer.processInput(); simulator.update(dt); renderer.render(simulator); renderer.swapBuffers(); - std::this_thread::sleep_for(std::chrono::milliseconds(16)); // Aim for roughly 60 FPS + std::this_thread::sleep_for(std::chrono::milliseconds(16)); } return 0;