Merge pull request #5 from Ben1980/develop

Develop

CMakeLists.txt

@@ -4,7 +4,7 @@ if(${CMAKE_VERSION} VERSION_LESS 3.12)
     cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
 endif()
 
-project(gravity VERSION 0.2.0
+project(gravity VERSION 0.4.0
     DESCRIPTION "N-Body-Problem project of https://thoughts-on-cpp.com"
     LANGUAGES CXX)
 

solver/CMakeLists.txt

@@ -1,5 +1,11 @@
 add_library(solver
-    src/solver.cpp include/particle.h include/particleBuilder.h src/particleBuilder.cpp include/types.h)
+        src/solver.cpp
+        include/particle.h
+        include/particleBuilder.h
+        src/particleBuilder.cpp
+        include/types.h
+        src/types.cpp
+        src/particle.cpp)
 
 target_include_directories(solver 
     PUBLIC 

solver/include/particle.h

@@ -2,16 +2,38 @@
 #define PARTICLE_H
 
 #include "types.h"
+#include <cstdlib>
 
-struct Particle {
+class Particle {
+public:
+    Particle();
+    Particle(double mass, const Vector2D &acceleration, const Vector2D &velocity, const Vector2D &position);
+
+    bool operator==(const Particle &rhs) const;
+    bool operator!=(const Particle &rhs) const;
+
+    double getMass() const;
+
+    const Vector2D &getAcceleration() const;
+
+    void setAcceleration(const Vector2D &acceleration);
+
+    const Vector2D &getVelocity() const;
+
+    void setVelocity(const Vector2D &velocity);
+
+    const Vector2D &getPosition() const;
+
+    void setPosition(const Vector2D &position);
+
+private:
+    static size_t IDCounter;
+
+    size_t id;
     double mass;
     Vector2D acceleration;
     Vector2D velocity;
     Vector2D position;
-
-    Particle() : mass(0) {}
-    Particle(double mass, const Vector2D &acceleration, const Vector2D &velocity, const Vector2D &position)
-        : mass(mass), acceleration(acceleration), velocity(velocity), position(position) {}
 };
 
 #endif

solver/include/particleBuilder.h

@@ -2,6 +2,7 @@
 #define PARTICLEBUILDER_H
 
 #include "types.h"
+#include <vector>
 
 class Particle;
 
@@ -14,6 +15,7 @@ class ParticleBuilder {
     ParticleBuilder & acceleration(const Vector2D &acceleration);
     ParticleBuilder & mass(double mass);
     Particle build() const;
+    std::vector<Particle> build(size_t numberOfParticles);
 
 private:
     double mMass;

solver/include/solver.h

@@ -2,12 +2,26 @@
 #define SOLVER_H
 
 #include <vector>
+#include "types.h"
 
 class Particle;
 
 class Solver {
 public:
-    std::vector<Particle> solve(const std::vector<Particle> &particles, double epsilon) const;
+    explicit Solver(double mEpsilon);
+
+    std::vector<Particle> solve(const std::vector<Particle> &particles) const;
+
+private:
+    std::vector<Particle> calculateAcceleration(const std::vector<Particle> &particles) const;
+    std::vector<Particle> calculateVelocity(const std::vector<Particle> &particles) const;
+    std::vector<Particle> calculatePosition(const std::vector<Particle> &particles) const;
+    static Particle AccumulateAcceleration(const std::vector<Particle> &particles, const Particle &particle);
+    static double CalculateEquivalentMass(const Particle &particleA, const Particle &particleB);
+
+    static const double G;
+    static const double EPSILON;
+    double mEpsilon;
 };
 
 #endif

solver/include/types.h

@@ -7,6 +7,30 @@ struct Vector2D {
 
     Vector2D() : x(0), y(0) {}
     Vector2D(double x, double y) : x(x), y(y) {}
+
+    bool operator==(const Vector2D &rhs) const;
+
+    bool operator!=(const Vector2D &rhs) const;
+
+    double length() const;
+
+    Vector2D& operator-=(const Vector2D& rhs);
+
+    Vector2D& operator+=(const Vector2D& rhs);
+
+    Vector2D& operator*=(const double& rhs);
+
+    Vector2D& operator/=(const double& rhs);
 };
 
+Vector2D operator-(const Vector2D &lhs, const Vector2D &rhs);
+
+Vector2D operator+(const Vector2D &lhs, const Vector2D &rhs);
+
+Vector2D operator*(const Vector2D &lhs, const double &rhs);
+
+Vector2D operator*(const double &lhs, const Vector2D &rhs);
+
+Vector2D operator/(const Vector2D &lhs, const double &rhs);
+
 #endif

solver/src/particle.cpp

@@ -0,0 +1,54 @@
+#include "particle.h"
+#include <random>
+#include <functional>
+#include <chrono>
+#include <cassert>
+
+size_t Particle::IDCounter = 0;
+
+Particle::Particle()
+    : mass(0) {
+    id = IDCounter++;
+}
+
+Particle::Particle(double mass, const Vector2D &acceleration, const Vector2D &velocity, const Vector2D &position)
+        : mass(mass), acceleration(acceleration), velocity(velocity), position(position) {
+    assert(mass > 0);
+    id = IDCounter++;
+}
+
+bool Particle::operator==(const Particle &rhs) const {
+    return id == rhs.id;
+}
+
+bool Particle::operator!=(const Particle &rhs) const {
+    return !(rhs == *this);
+}
+
+const Vector2D &Particle::getAcceleration() const {
+    return acceleration;
+}
+
+void Particle::setAcceleration(const Vector2D &acceleration) {
+    Particle::acceleration = acceleration;
+}
+
+const Vector2D &Particle::getVelocity() const {
+    return velocity;
+}
+
+void Particle::setVelocity(const Vector2D &velocity) {
+    Particle::velocity = velocity;
+}
+
+const Vector2D &Particle::getPosition() const {
+    return position;
+}
+
+void Particle::setPosition(const Vector2D &position) {
+    Particle::position = position;
+}
+
+double Particle::getMass() const {
+    return mass;
+}

solver/src/particleBuilder.cpp

@@ -1,5 +1,6 @@
 #include "particleBuilder.h"
 #include "particle.h"
+#include <random>
 
 ParticleBuilder & ParticleBuilder::position(const Vector2D &position)
 {
@@ -29,3 +30,21 @@ Particle ParticleBuilder::build() const
 {
     return {mMass, mAcceleration, mVelocity, mPosition};
 }
+
+std::vector<Particle> ParticleBuilder::build(size_t numberOfParticles)
+{
+    std::vector<Particle> particle(numberOfParticles);
+
+    std::mt19937 mt(std::random_device{}());
+    std::uniform_real_distribution real_dist(1.0, static_cast<double>(numberOfParticles));
+
+    for(Particle &p : particle) {
+        p = mass(real_dist(mt))
+                .acceleration({ real_dist(mt), real_dist(mt) })
+                .velocity({ real_dist(mt), real_dist(mt) })
+                .position({ real_dist(mt), real_dist(mt) })
+                .build();
+    }
+
+    return particle;
+}

solver/src/solver.cpp

@@ -1,7 +1,85 @@
 #include "solver.h"
 #include "particle.h"
+#include <cmath>
+#include <algorithm>
+#include <solver.h>
+#include <cassert>
 
-std::vector<Particle> Solver::solve(const std::vector<Particle> &particles, double epsilon) const
-{
-    return std::vector<Particle>(2) ;
-}
+
+const double Solver::G = 6.67408e-11;
+const double Solver::EPSILON = 1e-3;
+
+Solver::Solver(double mEpsilon) : mEpsilon(mEpsilon) {}
+
+std::vector<Particle> Solver::solve(const std::vector<Particle> &particles) const {
+    std::vector<Particle> solution = calculateAcceleration(particles);
+    solution = calculateVelocity(solution);
+    solution = calculatePosition(solution);
+
+    return solution;
+}
+
+std::vector<Particle> Solver::calculateAcceleration(const std::vector<Particle> &particles) const {
+    std::vector<Particle> solution(particles.size());
+
+    std::transform(begin(particles), end(particles), begin(solution), [&particles](const Particle &particle) {
+        return AccumulateAcceleration(particles, particle);
+    });
+
+    return solution;
+}
+
+std::vector<Particle> Solver::calculateVelocity(const std::vector<Particle> &particles) const {
+    std::vector<Particle> solution(particles.size());
+
+    std::transform(begin(particles), end(particles), begin(solution), [epsilon = mEpsilon](Particle particle) {
+        const Vector2D v0 = particle.getVelocity();
+        particle.setVelocity(v0 + particle.getAcceleration()*epsilon);
+
+        return particle;
+    });
+
+    return solution;
+}
+
+std::vector<Particle> Solver::calculatePosition(const std::vector<Particle> &particles) const {
+    std::vector<Particle> solution(particles.size());
+
+    std::transform(begin(particles), end(particles), begin(solution), [epsilon = mEpsilon](Particle particle) {
+        const Vector2D v = particle.getVelocity();
+        const Vector2D r0 = particle.getPosition();
+        particle.setPosition(r0 + v*epsilon + particle.getAcceleration()*epsilon*epsilon);
+
+        return particle;
+    });
+
+    return solution;
+}
+
+Particle Solver::AccumulateAcceleration(const std::vector<Particle> &particles, const Particle &particle) {
+    Particle particleA = particle;
+    const double e3 = EPSILON*EPSILON*EPSILON;
+
+    std::for_each(begin(particles), end(particles), [&particleA, e3](const Particle &particleB) {
+        if(particleA != particleB) {
+            const double M = CalculateEquivalentMass(particleA, particleB);
+            const Vector2D r = particleB.getPosition() - particleA.getPosition();
+            const double rLength = r.length();
+            const double r3 = fabs(rLength*rLength*rLength);
+
+            const Vector2D a0 = particleA.getAcceleration();
+            particleA.setAcceleration(a0 + G*M*r/(r3 + e3));
+        }
+    });
+
+    return particleA;
+}
+
+double Solver::CalculateEquivalentMass(const Particle &particleA, const Particle &particleB) {
+    const double massA = particleA.getMass();
+    assert(massA > 0);
+
+    const double massB = particleB.getMass();
+
+    return massA*massB/massA;
+}

solver/src/types.cpp

@@ -0,0 +1,67 @@
+#include "types.h"
+
+#include <cmath>
+#include <types.h>
+#include <cassert>
+#include <limits>
+
+bool Vector2D::operator==(const Vector2D &rhs) const {
+    auto equalZero = std::numeric_limits<double>::min();
+
+    return fabs(x - rhs.x) <= equalZero &&
+           fabs(y - rhs.y) <= equalZero;
+}
+
+bool Vector2D::operator!=(const Vector2D &rhs) const {
+    return !(rhs == *this);
+}
+
+double Vector2D::length() const {
+    return sqrt(x*x + y*y);
+}
+
+Vector2D& Vector2D::operator-=(const Vector2D& rhs) {
+    *this = *this - rhs;
+    return *this;
+}
+
+Vector2D& Vector2D::operator*=(const double& rhs) {
+    *this = *this * rhs;
+    return *this;
+}
+
+Vector2D& Vector2D::operator/=(const double& rhs) {
+    *this = *this / rhs;
+    return *this;
+}
+
+Vector2D &Vector2D::operator+=(const Vector2D &rhs) {
+    *this = *this + rhs;
+    return *this;
+}
+
+Vector2D operator-(const Vector2D &lhs, const Vector2D &rhs) {
+    return { lhs.x - rhs.x,
+             lhs.y - rhs.y };
+}
+
+Vector2D operator+(const Vector2D &lhs, const Vector2D &rhs) {
+    return { lhs.x + rhs.x,
+             lhs.y + rhs.y };
+}
+
+Vector2D operator*(const Vector2D &lhs, const double &rhs) {
+    return { lhs.x * rhs,
+             lhs.y * rhs };
+}
+
+Vector2D operator*(const double &lhs, const Vector2D &rhs) {
+    return rhs * lhs;
+}
+
+Vector2D operator/(const Vector2D &lhs, const double &rhs) {
+    assert(fabs(rhs) > 0);
+
+    return { lhs.x / rhs,
+             lhs.y / rhs };
+}