diff --git a/CelestialBody.cpp b/CelestialBody.cpp index c6ab7f0..66f0f35 100644 --- a/CelestialBody.cpp +++ b/CelestialBody.cpp @@ -12,24 +12,28 @@ std::string vec3_to_string(const glm::vec3& v) { return ss.str(); } -CelestialBody::CelestialBody(double mass, const glm::dvec3& position, const glm::dvec3& velocity) +CelestialBody::CelestialBody(double mass, const glm::dvec3& position, const glm::dvec3& velocity, double radius) : mass(mass), position(position), velocity(velocity), acceleration(0.0f) {} 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); - } + // Runge-Kutta 4th order method + glm::dvec3 k1v = acceleration * dt; + glm::dvec3 k1r = velocity * dt; - velocity += acceleration * dt; - position += velocity * dt; + glm::dvec3 k2v = acceleration * dt; + glm::dvec3 k2r = (velocity + k1v * 0.5) * dt; - if (glm::any(glm::isnan(position)) || glm::any(glm::isinf(position))) { - std::cout << "Warning: Invalid position detected: " << vec3_to_string(position) << std::endl; - position = glm::vec3(0.0f); - } + glm::dvec3 k3v = acceleration * dt; + glm::dvec3 k3r = (velocity + k2v * 0.5) * dt; + + glm::dvec3 k4v = acceleration * dt; + glm::dvec3 k4r = (velocity + k3v) * dt; + + velocity += (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0; + position += (k1r + 2.0 * k2r + 2.0 * k3r + k4r) / 6.0; + + acceleration = glm::dvec3(0.0); - acceleration = glm::vec3(0.0f); } void CelestialBody::applyForce(const glm::dvec3& force) { diff --git a/CelestialBody.h b/CelestialBody.h index 1c22a60..a4e7d20 100644 --- a/CelestialBody.h +++ b/CelestialBody.h @@ -11,7 +11,7 @@ class CelestialBody { public: - CelestialBody(double mass, const glm::dvec3& position, const glm::dvec3& velocity); + CelestialBody(double mass, const glm::dvec3& position, const glm::dvec3& velocity, double radius); void update(double dt); void applyForce(const glm::dvec3& force); @@ -21,6 +21,7 @@ public: [[nodiscard]] glm::dvec3 getVelocity() const { return velocity; } void addToTrajectory(const glm::dvec3& position); const std::vector& getTrajectory() const { return trajectory; } + double getRadius() const { return radius; } private: double mass; @@ -29,5 +30,6 @@ private: glm::dvec3 acceleration; std::vector trajectory; static const size_t MAX_TRAJECTORY_POINTS = 1000; + double radius; }; #endif //GRAVITY_CELESTIALBODY_H diff --git a/Renderer.cpp b/Renderer.cpp index d6ea9bf..fa5d9eb 100644 --- a/Renderer.cpp +++ b/Renderer.cpp @@ -2,12 +2,9 @@ // Created by Quinta on 7/12/2024. // #include "Renderer.h" -#include -#include #include #include #include -#include Renderer::Renderer(int width, int height) : cameraPos(3e11f, 2e11f, 3e11f), @@ -65,7 +62,7 @@ void Renderer::render(const Simulator& simulator) { glMatrixMode(GL_PROJECTION); glLoadIdentity(); - gluPerspective(45.0, 1600.0 / 1200.0, 1e9, 1e13); + gluPerspective(45.0, 1600.0 / 1200.0, 1e8, 1e14); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); @@ -90,13 +87,9 @@ void Renderer::render(const Simulator& simulator) { 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]; glm::dvec3 pos = body.getPosition(); - //std::cout << "Body " << i << " position: " << pos.x << ", " << pos.y << ", " << pos.z << std::endl; } // Calculate the log range @@ -117,19 +110,6 @@ void Renderer::render(const Simulator& simulator) { glm::dvec3 pos = body.getPosition(); glm::vec3 renderPos(static_cast(pos.x), static_cast(pos.y), static_cast(pos.z)); -// std::cout << "Rendering body " << i << " ("; -// switch(i) { -// case 0: std::cout << "Sun"; break; -// case 1: std::cout << "Mercury"; break; -// case 2: std::cout << "Venus"; break; -// case 3: std::cout << "Earth"; break; -// case 4: std::cout << "Mars"; break; -// default: std::cout << "Unknown"; break; -// } -// std::cout << ") at position (" -// << pos.x << ", " << pos.y << ", " << pos.z -// << ") with scale " << scaleFactor << std::endl; - // Set color based on body index switch(i) { case 0: glColor3f(1.0f, 1.0f, 0.0f); break; // Sun: Yellow @@ -137,6 +117,12 @@ void Renderer::render(const Simulator& simulator) { case 2: glColor3f(0.9f, 0.7f, 0.4f); break; // Venus: Light Orange case 3: glColor3f(0.0f, 0.5f, 1.0f); break; // Earth: Blue case 4: glColor3f(1.0f, 0.0f, 0.0f); break; // Mars: Red + case 5: glColor3f(0.8f, 0.6f, 0.2f); break; // Jupiter: Light Brown + case 6: glColor3f(0.9f, 0.9f, 0.7f); break; // Saturn: Light Yellow + case 7: glColor3f(0.0f, 0.5f, 0.5f); break; // Uranus: Cyan + case 8: glColor3f(0.0f, 0.0f, 1.0f); break; // Neptune: Dark Blue + case 9: glColor3f(0.5f, 0.5f, 0.5f); break; // Pluto: Gray + default: glColor3f(1.0f, 1.0f, 1.0f); break; // White for any additional bodies } @@ -271,8 +257,8 @@ float Renderer::calculateGravityFieldStrength(const glm::vec3& point, const std: } void Renderer::drawGrid(const Simulator& simulator) { - const float gridSize = 5e11f; - const int gridLines = 20; + const float gridSize = 5e13f; + const int gridLines = 80; const float lineSpacing = gridSize / gridLines; glBegin(GL_LINES); @@ -311,7 +297,7 @@ void Renderer::drawTrajectories(const std::vector& bodies) { } void Renderer::processInput() { - float cameraSpeed = this->cameraSpeed; + float cameraSpeed = this->cameraSpeed * 1e1f; glm::vec3 front(cameraFront.x, 0, cameraFront.z); front = glm::normalize(front); diff --git a/Simulator.cpp b/Simulator.cpp index f439ff7..e09ffbe 100644 --- a/Simulator.cpp +++ b/Simulator.cpp @@ -35,6 +35,9 @@ void Simulator::update(double dt) { bodies[i].update(dt); bodies[i].addToTrajectory(bodies[i].getPosition()); } + + // Check for collisions + checkCollisions(); } glm::dvec3 Simulator::calculateGravitationalForce(const CelestialBody& body1, const CelestialBody& body2) { @@ -52,10 +55,49 @@ glm::dvec3 Simulator::calculateGravitationalForce(const CelestialBody& body1, co 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::dvec3(0.0); } return glm::normalize(direction) * forceMagnitude; +} + +void Simulator::handleCollision(CelestialBody& body1, CelestialBody& body2) { + double totalMass = body1.getMass() + body2.getMass(); + + // Calculate center of mass position + glm::dvec3 newPosition = (body1.getPosition() * body1.getMass() + body2.getPosition() * body2.getMass()) / totalMass; + + // Calculate new velocity (momentum conservation) + glm::dvec3 newVelocity = (body1.getVelocity() * body1.getMass() + body2.getVelocity() * body2.getMass()) / totalMass; + + // Calculate new radius (assuming constant density) + double newRadius = std::pow(std::pow(body1.getRadius(), 3) + std::pow(body2.getRadius(), 3), 1.0/3.0); + + // Create new body + CelestialBody newBody(totalMass, newPosition, newVelocity, newRadius); + + // Replace body1 with the new body + body1 = newBody; + + // Remove body2 + auto it = std::find_if(bodies.begin(), bodies.end(), [&body2](const CelestialBody& b) { + return &b == &body2; + }); + if (it != bodies.end()) { + bodies.erase(it); + } +} + +void Simulator::checkCollisions() { + for (size_t i = 0; i < bodies.size(); ++i) { + for (size_t j = i + 1; j < bodies.size(); ++j) { + CelestialBody& body1 = bodies[i]; + CelestialBody& body2 = bodies[j]; + glm::dvec3 distanceVec = body1.getPosition() - body2.getPosition(); + double distance = glm::length(distanceVec); + if (distance < (body1.getRadius() + body2.getRadius())) { + handleCollision(body1, body2); + } + } + } } \ No newline at end of file diff --git a/Simulator.h b/Simulator.h index 8d40590..0852f55 100644 --- a/Simulator.h +++ b/Simulator.h @@ -21,5 +21,7 @@ public: private: std::vector bodies; const float G = 6.67430e-11f; // Gravitational constant + void checkCollisions(); + void handleCollision(CelestialBody& body1, CelestialBody& body2); }; #endif //GRAVITY_SIMULATOR_H diff --git a/main.cpp b/main.cpp index 3eb79ab..788e1b0 100644 --- a/main.cpp +++ b/main.cpp @@ -6,24 +6,58 @@ #include #include +glm::dvec3 calculateOrbitalVelocity(double centralMass, double distance) { + const double G = 6.67430e-11; + double speed = std::sqrt(G * centralMass / distance); + return glm::dvec3(0, speed, 0); // Assuming orbit in the XZ plane +} + int main() { Simulator simulator; Renderer renderer(1600, 1200); -// Sun (at the center) - simulator.addBody(CelestialBody(1.989e30f, glm::dvec3(0, 0, 0), glm::dvec3(0, 0, 0))); + double sunMass = 1.989e30; -// Mercury - simulator.addBody(CelestialBody(3.285e23f, glm::dvec3(57.9e9f, 0, 0), glm::dvec3(0, 47.36e3f, 0))); + // Sun (at the center) + simulator.addBody(CelestialBody(sunMass, glm::dvec3(0, 0, 0), glm::dvec3(0, 0, 0), 6.96e8)); -// Venus - simulator.addBody(CelestialBody(4.867e24f, glm::dvec3(108.2e9f, 0, 0), glm::dvec3(0, 35.02e3f, 0))); + // Mercury + double mercuryDist = 57.9e9; + simulator.addBody(CelestialBody(3.285e23, glm::dvec3(mercuryDist, 0, 0), calculateOrbitalVelocity(sunMass, mercuryDist), 2.44e6)); + + // Venus + double venusDist = 108.2e9; + simulator.addBody(CelestialBody(4.867e24, glm::dvec3(venusDist, 0, 0), calculateOrbitalVelocity(sunMass, venusDist), 6.05e6)); + + // Earth + double earthDist = 149.6e9; + simulator.addBody(CelestialBody(5.972e24, glm::dvec3(earthDist, 0, 0), calculateOrbitalVelocity(sunMass, earthDist), 6.37e6)); + + // Mars + double marsDist = 227.9e9; + simulator.addBody(CelestialBody(6.39e23, glm::dvec3(marsDist, 0, 0), calculateOrbitalVelocity(sunMass, marsDist), 3.39e6)); + + // Jupiter + double jupiterDist = 778.5e9; + simulator.addBody(CelestialBody(1.898e27, glm::dvec3(jupiterDist, 0, 0), calculateOrbitalVelocity(sunMass, jupiterDist), 69.91e6)); + + // Saturn + double saturnDist = 1.429e12; + simulator.addBody(CelestialBody(5.683e26, glm::dvec3(saturnDist, 0, 0), calculateOrbitalVelocity(sunMass, saturnDist), 58.23e6)); + + // Uranus + double uranusDist = 2.871e12; + simulator.addBody(CelestialBody(8.681e25, glm::dvec3(uranusDist, 0, 0), calculateOrbitalVelocity(sunMass, uranusDist ), 25.36e6)); + + // Neptune + double neptuneDist = 4.495e12; + simulator.addBody(CelestialBody(1.024e26, glm::dvec3(neptuneDist, 0, 0), calculateOrbitalVelocity(sunMass, neptuneDist), 24.62e6)); + + // Pluto + double plutoDist = 5.906e12; + simulator.addBody(CelestialBody(1.309e22, glm::dvec3(plutoDist, 0, 0), calculateOrbitalVelocity(sunMass, plutoDist), 1.18e6)); -// Earth - simulator.addBody(CelestialBody(5.972e24f, glm::dvec3(149.6e9f, 0, 0), glm::dvec3(0, 29.78e3f, 0))); -// Mars - 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