Merge pull request 'SImplify acceleration addition.' (#1) from tfaour-patch-1 into main

Reviewed-on: #1

CMakeLists.txt

@@ -1,12 +1,17 @@
 cmake_minimum_required(VERSION 3.10)
 project(orbital_simulator)
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 23)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
+# Enable testing
+enable_testing()
+include(CTest)
+
 # Find required packages
 find_package(Boost REQUIRED COMPONENTS program_options)
 find_package(nlohmann_json REQUIRED)
+find_package(GTest REQUIRED)
 
 # Add source files
 set(SOURCES
@@ -16,6 +21,13 @@ set(SOURCES
     src/simulator.cpp
 )
 
+# Add source files for tests (excluding main.cpp)
+set(TEST_SOURCES
+    src/body.cpp
+    src/calc.cpp
+    src/simulator.cpp
+)
+
 # Add header files
 set(HEADERS
     src/body.hpp
@@ -24,22 +36,46 @@ set(HEADERS
     src/units.hpp
 )
 
-# Create executable
+# Create main executable
 add_executable(orbital_simulator ${SOURCES} ${HEADERS})
 
-# Link libraries
+# Create test executable
+add_executable(orbital_simulator_tests
+    tests/body_test.cpp
+    tests/calc_test.cpp
+    tests/simulator_test.cpp
+    ${TEST_SOURCES}  # Use test sources instead of all sources
+)
+
+# Link libraries for main executable
 target_link_libraries(orbital_simulator
     PRIVATE
     Boost::program_options
     nlohmann_json::nlohmann_json
 )
 
+# Link libraries for test executable
+target_link_libraries(orbital_simulator_tests
+    PRIVATE
+    GTest::GTest
+    GTest::Main
+)
+
 # Include directories
 target_include_directories(orbital_simulator
     PRIVATE
     ${CMAKE_SOURCE_DIR}/src
 )
 
+target_include_directories(orbital_simulator_tests
+    PRIVATE
+    ${CMAKE_SOURCE_DIR}/src
+    ${CMAKE_SOURCE_DIR}/tests
+)
+
+# Add tests to CTest
+add_test(NAME orbital_simulator_tests COMMAND orbital_simulator_tests)
+
 # Optional: Enable ncurses for terminal plotting
 option(ENABLE_NCURSES "Enable terminal plotting with ncurses" OFF)
 if(ENABLE_NCURSES)

src/body.cpp

@@ -1,13 +1,20 @@
 #include "body.hpp"
 #include "calc.hpp"
 #include <cmath>
-#include <boost/multiprecision/cpp_dec_float.hpp>
+#include <ranges>
+
 
 Body::Body(const Position& X, const Velocity& V, const Mass& m, const std::string& name)
     : X(X), V(V), m(m), name(name) {
     A = Acceleration{Decimal(0), Decimal(0), Decimal(0)};
 }
 
+void Body::addAcceleration(const Acceleration& new_A) {
+    for (const auto& [new_a, cur_a]: std::views::zip(new_A, this->A)){
+        cur_a += new_a;
+    }
+}
+
 std::tuple<Position, Velocity, Mass> Body::save() const {
     return std::make_tuple(X, V, m);
 }

src/body.hpp

@@ -30,6 +30,8 @@ public:
     // Setters
     void setAcceleration(const Acceleration& new_A) { A = new_A; }
     
+    void addAcceleration(const Acceleration& new_A);
+    
     // String representation
     std::string toString() const;
 

src/simulator.cpp

@@ -88,8 +88,8 @@ void Simulator::calculate_forces() {
             
             // Calculate force magnitude using Newton's law of gravitation
             // F = G * m1 * m2 / r^2 BUT G = 1, and we'll multiply by the opposite mass later
-            // for the acceleration
-            Decimal force_magnitude = 1 / (dist * dist);
+            // for the acceleration. Use r^3 to avoid normalizing vec when multiplying later
+            Decimal force_magnitude = 1 / (dist * dist * dist); 
             
             // Calculate acceleration for both bodies
             Decimal acc_magnitude_i = force_magnitude * bodies[j].getMass();
@@ -98,19 +98,13 @@ void Simulator::calculate_forces() {
             // Convert to vector form
             Acceleration acc_i, acc_j;
             for (int k = 0; k < 3; ++k) {
-                acc_i[k] = -vec[k] * acc_magnitude_i / dist;  
-                acc_j[k] = vec[k] * acc_magnitude_j / dist;   
+                acc_i[k] = -vec[k] * acc_magnitude_i;  
+                acc_j[k] = vec[k] * acc_magnitude_j;   
             }
             
-            // Add to current accelerations
-            Acceleration current_acc_i = bodies[i].getAcceleration();
-            Acceleration current_acc_j = bodies[j].getAcceleration();
-            for (int k = 0; k < 3; ++k) {
-                current_acc_i[k] += acc_i[k];
-                current_acc_j[k] += acc_j[k];
-            }
-            bodies[i].setAcceleration(current_acc_i);
-            bodies[j].setAcceleration(current_acc_j);
+            // Add this component to the two bodies
+            bodies[i].addAcceleration(acc_i);
+            bodies[j].addAcceleration(acc_j);
         }
     }
 }

src/simulator.hpp

@@ -22,6 +22,13 @@ public:
     // Run simulation
     void run(int steps);
 
+    // Getters
+    const std::vector<Body>& getBodies() const { return bodies; }
+    Decimal getStepSize() const { return step_size; }
+    int getStepsPerSave() const { return steps_per_save; }
+    const std::filesystem::path& getOutputFile() const { return output_file; }
+    int getCurrentStep() const { return current_step; }
+
 private:
     void calculate_forces();
     void move_bodies();

src/units.hpp

@@ -41,4 +41,6 @@ inline Decimal real_mass(Decimal mass) { return mass * m_0; }
 inline Decimal norm_time(Decimal time) { return time / t_0; }
 inline Decimal real_time(Decimal time) { return time * t_0; }
 inline Decimal norm_vel(Decimal vel) { return vel / (r_0/t_0); }
-inline Decimal real_vel(Decimal vel) { return vel * (r_0/t_0); } 
+inline Decimal real_vel(Decimal vel) { return vel * (r_0/t_0); }
+inline Decimal norm_acc(Decimal acc) { return acc / (r_0/(t_0*t_0)); }
+inline Decimal real_acc(Decimal acc) { return acc * (r_0/(t_0*t_0)); } 

tests/body_test.cpp

@@ -0,0 +1,111 @@
+#include <gtest/gtest.h>
+#include "body.hpp"
+#include <cmath>
+#include <memory>
+
+class BodyTest : public testing::Test {
+protected:
+    BodyTest() {
+        test_body = std::make_unique<Body>(
+            Position{Decimal(0), Decimal(0), Decimal(0)},
+            Velocity{Decimal(0), Decimal(0), Decimal(0)},
+            Mass(1.0),
+            "test_body"
+        );
+    }
+
+    std::unique_ptr<Body> test_body;
+};
+
+TEST_F(BodyTest, AddAcceleration) {
+    Acceleration new_A{Decimal(1.0), Decimal(2.0), Decimal(3.0)};
+    test_body->addAcceleration(new_A);
+
+    const auto& A = test_body->getAcceleration();
+    EXPECT_DOUBLE_EQ(A[0], 1.0);
+    EXPECT_DOUBLE_EQ(A[1], 2.0);
+    EXPECT_DOUBLE_EQ(A[2], 3.0);
+
+    // Add another acceleration
+    Acceleration another_A{Decimal(2.0), Decimal(3.0), Decimal(4.0)};
+    test_body->addAcceleration(another_A);
+
+    EXPECT_DOUBLE_EQ(A[0], 3.0);
+    EXPECT_DOUBLE_EQ(A[1], 5.0);
+    EXPECT_DOUBLE_EQ(A[2], 7.0);
+}
+
+TEST_F(BodyTest, Step) {
+    // Set initial velocity and acceleration
+    test_body->setAcceleration(Acceleration{Decimal(1.0), Decimal(2.0), Decimal(3.0)});
+    
+    // Take a step
+    test_body->step(1.0);
+
+    // Check new position
+    const auto& X = test_body->getPosition();
+    EXPECT_DOUBLE_EQ(X[0], 0.0);
+    EXPECT_DOUBLE_EQ(X[1], 0.0);
+    EXPECT_DOUBLE_EQ(X[2], 0.0);
+
+    // Check new velocity
+    const auto& V = test_body->getVelocity();
+    EXPECT_DOUBLE_EQ(V[0], 1.0);
+    EXPECT_DOUBLE_EQ(V[1], 2.0);
+    EXPECT_DOUBLE_EQ(V[2], 3.0);
+
+    // Check acceleration is reset
+    const auto& A = test_body->getAcceleration();
+    EXPECT_DOUBLE_EQ(A[0], 0.0);
+    EXPECT_DOUBLE_EQ(A[1], 0.0);
+    EXPECT_DOUBLE_EQ(A[2], 0.0);
+
+}
+
+TEST_F(BodyTest, Energy) {
+    // Set position and velocity through save/load
+    auto state = std::make_tuple(
+        Position{Decimal(1.0), Decimal(0.0), Decimal(0.0)},
+        Velocity{Decimal(1.0), Decimal(0.0), Decimal(0.0)},
+        Mass(1.0)
+    );
+    test_body = std::make_unique<Body>(Body::load(state));
+
+    // Calculate expected values
+    Decimal expected_ke = Decimal(0.5); // 0.5 * m * v^2
+    Decimal expected_pe = Decimal(-1.0); // -m/r
+    Decimal expected_total = expected_ke + expected_pe;
+
+    EXPECT_DOUBLE_EQ(test_body->ke(), expected_ke);
+    EXPECT_DOUBLE_EQ(test_body->pe(), expected_pe);
+    EXPECT_DOUBLE_EQ(test_body->E(), expected_total);
+}
+
+TEST_F(BodyTest, SaveAndLoad) {
+    // Set some values through save/load
+    auto state = std::make_tuple(
+        Position{Decimal(1.0), Decimal(2.0), Decimal(3.0)},
+        Velocity{Decimal(4.0), Decimal(5.0), Decimal(6.0)},
+        Mass(7.0)
+    );
+    test_body = std::make_unique<Body>(Body::load(state));
+
+    // Save state
+    auto saved_state = test_body->save();
+
+    // Create new body from saved state
+    Body loaded_body = Body::load(saved_state);
+
+    // Verify loaded values
+    const auto& X = loaded_body.getPosition();
+    const auto& V = loaded_body.getVelocity();
+    const auto& m = loaded_body.getMass();
+
+    EXPECT_DOUBLE_EQ(X[0], 1.0);
+    EXPECT_DOUBLE_EQ(X[1], 2.0);
+    EXPECT_DOUBLE_EQ(X[2], 3.0);
+    EXPECT_DOUBLE_EQ(V[0], 4.0);
+    EXPECT_DOUBLE_EQ(V[1], 5.0);
+    EXPECT_DOUBLE_EQ(V[2], 6.0);
+    EXPECT_DOUBLE_EQ(m, 7.0);
+} 

tests/calc_test.cpp

@@ -0,0 +1,10 @@
+#include <gtest/gtest.h>
+#include "calc.hpp"
+#include <cmath>
+
+TEST(CalcTest, FormatSigFigs) {
+    // Test positive numbers
+    EXPECT_EQ(format_sig_figs(123.456, 3), "1.23e+02");
+    EXPECT_EQ(format_sig_figs(123.456, 4), "1.235e+02");
+    EXPECT_EQ(format_sig_figs(123.456, 5), "1.2346e+02");
+}

tests/simulator_test.cpp

@@ -0,0 +1,54 @@
+#include <gtest/gtest.h>
+#include "simulator.hpp"
+#include <vector>
+#include <filesystem>
+#include <memory>
+
+class SimulatorTest : public testing::Test {
+protected:
+    SimulatorTest() : 
+        bodies({
+            Body(
+                Position{Decimal(0), Decimal(0), Decimal(0)},
+                Velocity{Decimal(0), Decimal(0), Decimal(0)},
+                Mass(1.0),
+                "body1"
+            ),
+            Body(
+                Position{Decimal(1), Decimal(0), Decimal(0)},
+                Velocity{Decimal(0), Decimal(1), Decimal(0)},
+                Mass(1.0),
+                "body2"
+            )
+        })
+    {
+        simulator = std::make_unique<Simulator>(
+            bodies,
+            Decimal(0.1),  // step_size
+            1,            // steps_per_save
+            std::filesystem::path("test_output.json"),  // output_file
+            0,            // current_step
+            true          // overwrite_output
+        );
+    }
+
+    std::vector<Body> bodies;
+    std::unique_ptr<Simulator> simulator;
+};
+
+TEST_F(SimulatorTest, Step) {
+    // Take a step
+    //TODO
+}
+
+TEST_F(SimulatorTest, TotalEnergy) {
+    //TODO
+    // Decimal initial_energy = simulator->total_energy();
+
+    // simulator->step(Decimal(0.1));
+
+    // Decimal new_energy = simulator->total_energy();
+
+    // EXPECT_NEAR(initial_energy, new_energy, 1e-10);
+}
+