Resolved that annoying gravity bug and added trajectories, a grid, keyboard navigation and logarithimic scale

.idea/vcs.xml

@@ -1,6 +1,6 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
   <component name="VcsDirectoryMappings">
-    <mapping directory="$PROJECT_DIR$" vcs="Git" />
+    <mapping directory="" vcs="Git" />
   </component>
 </project>

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

CelestialBody.h

@@ -7,22 +7,27 @@
 #pragma once
 #include <glm/glm.hpp>
 #include <string>
+#include <vector>
 
 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<glm::dvec3>& 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<glm::dvec3> trajectory;
+    static const size_t MAX_TRAJECTORY_POINTS = 1000;
 };
 #endif //GRAVITY_CELESTIALBODY_H

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.
+

Renderer.cpp

@@ -9,7 +9,18 @@
 #include <stdexcept>
 #include <iostream>
 
-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<Renderer*>(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<double>::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<float>(normalizedLogMass) * (maxScale - minScale);
+
+        glm::dvec3 pos = body.getPosition();
+        glm::vec3 renderPos(static_cast<float>(pos.x), static_cast<float>(pos.y), static_cast<float>(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<CelestialBody>& 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<float>(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<CelestialBody>& 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<float>(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);
 }

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<CelestialBody>& bodies);
+    void drawGravityField(const Simulator& simulator);
+    void drawTrajectories(const std::vector<CelestialBody>& 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

Simulator.cpp

@@ -4,6 +4,7 @@
 #include "Simulator.h"
 #include <glm/glm.hpp>
 #include <iostream>
+#include <algorithm>
 
 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;

Simulator.h

@@ -13,13 +13,13 @@ public:
     Simulator();
 
     void addBody(const CelestialBody& body);
-    void update(float dt);
+    void update(double dt);
     const std::vector<CelestialBody>& getBodies() const { return bodies; }
+    glm::dvec3 calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2);
+
 
 private:
     std::vector<CelestialBody> bodies;
     const float G = 6.67430e-11f; // Gravitational constant
-
-    glm::vec3 calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2);
 };
 #endif //GRAVITY_SIMULATOR_H

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;