Merge branch 'main' into tfaour-patch-1

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)

animate.py

@@ -50,19 +50,45 @@ def read_output_file(filename):
                 
             name = data['name']
             pos = data['position']
+            
+            # Store position
             positions[name].append((pos[0], pos[1]))
             frame += 1
             times.append(frame)
     
     return positions, times
 
-def create_animation(positions, times, output_file=None):
+def center_positions(positions, center_body):
+    """Center all positions relative to the specified body."""
+    if center_body not in positions:
+        print(f"Warning: Center body '{center_body}' not found in simulation")
+        return positions
+        
+    centered_positions = defaultdict(list)
+    for frame in range(len(next(iter(positions.values())))):
+        # Get center body position for this frame
+        center_pos = positions[center_body][frame]
+        
+        # Shift all positions relative to center body
+        for name, pos_list in positions.items():
+            if frame < len(pos_list):
+                pos = pos_list[frame]
+                centered_pos = (pos[0] - center_pos[0], pos[1] - center_pos[1])
+                centered_positions[name].append(centered_pos)
+    
+    return centered_positions
+
+def create_animation(positions, times, output_file=None, center_body=None):
     """Create an animation of the bodies' orbits."""
     # Check if we have any data
     if not positions or not times:
         print("Error: No valid data found in the input file")
         return
         
+    # Center positions if requested
+    if center_body:
+        positions = center_positions(positions, center_body)
+        
     # Set up the figure and axis
     fig, ax = plt.subplots(figsize=(10, 10))
     
@@ -77,7 +103,10 @@ def create_animation(positions, times, output_file=None):
     # Set up the plot
     ax.set_xlabel('X (m)')
     ax.set_ylabel('Y (m)')
-    ax.set_title('Orbital Simulation')
+    title = 'Orbital Simulation'
+    if center_body:
+        title += f' (centered on {center_body})'
+    ax.set_title(title)
     ax.legend()
     
     # Find the bounds of the plot
@@ -126,10 +155,11 @@ def main():
     parser = argparse.ArgumentParser(description='Animate orbital simulation output')
     parser.add_argument('input_file', help='Input file from simulation')
     parser.add_argument('--output', '-o', help='Output video file (optional)')
+    parser.add_argument('--center', '-c', help='Center the animation on this body')
     args = parser.parse_args()
     
     positions, times = read_output_file(args.input_file)
-    create_animation(positions, times, args.output)
+    create_animation(positions, times, args.output, args.center)
 
 if __name__ == '__main__':
     main() 

config/planets.json

@@ -1,16 +1,28 @@
 {
-    "planets": [
+    "bodies": [
+        {
+            "name": "Mercury",
+            "mass": 3.30104e23,
+            "position": [4.6000e10, 0, 0],
+            "velocity": [0, 58970, 0]
+        },
+        {
+            "name": "Venus",
+            "mass": 4.867e24,
+            "position": [1.08941e11, 0, 0],
+            "velocity": [0, 34780, 0]
+        },
         {
             "name": "Earth",
             "mass": 5.972e24,
-            "position": [149597870700, 0, 0],
-            "velocity": [0, 29780, 0]
+            "position": [1.47095e11, 0, 0],
+            "velocity": [0, 29290, 0]
         },
         {
             "name": "Moon",
             "mass": 7.34767309e22,
             "position": [149982270700, 0, 0],
-            "velocity": [0, 30802, 0]
+            "velocity": [0, 30822, 0]
         },
         {
             "name": "Sun",

src/body.cpp

@@ -1,13 +1,25 @@
 #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;
+    }
+}
+
+void Body::subAcceleration(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);
 }
@@ -59,8 +71,8 @@ std::string Body::speed() const {
 std::string Body::toString() const {
     std::string pos_str, vel_str;
     for (int i = 0; i < 3; ++i) {
-        pos_str += format_sig_figs(real_pos(X[i]), 3) + "m ";
-        vel_str += format_sig_figs(real_vel(V[i]), 3) + "m/s ";
+        pos_str += format_sig_figs(real_pos(X[i]), 10) + "m ";
+        vel_str += format_sig_figs(real_vel(V[i]), 10) + "m/s ";
     }
     return name + ": X = " + pos_str + ", V = " + vel_str;
 } 

src/body.hpp

@@ -3,10 +3,12 @@
 #include "units.hpp"
 #include <string>
 #include <tuple>
+#include <memory>
 
 class Body {
 public:
-    Body(const Position& X, const Velocity& V, const Mass& m, const std::string& name = "");
+    Body(const Position& X, const Velocity& V, const Mass& m, 
+         const std::string& name = "");
     
     // Save and load state
     std::tuple<Position, Velocity, Mass> save() const;
@@ -21,8 +23,8 @@ public:
     std::string speed() const;  // Speed as formatted string
     
     // Getters
-    const Position& getPosition() const { return X; }
-    const Velocity& getVelocity() const { return V; }
+    const Position& getPosition() const { return X; } 
+    const Velocity& getVelocity() const { return V; } 
     const Acceleration& getAcceleration() const { return A; }
     const Mass& getMass() const { return m; }
     const std::string& getName() const { return name; }
@@ -30,12 +32,14 @@ public:
     // Setters
     void setAcceleration(const Acceleration& new_A) { A = new_A; }
     
+    void addAcceleration(const Acceleration& new_A);
+    void subAcceleration(const Acceleration& new_A);
     // String representation
     std::string toString() const;
 
 private:
-    Position X;
-    Velocity V;
+    Position X; 
+    Velocity V;  
     Acceleration A;
     Mass m;
     std::string name;

src/main.cpp

@@ -52,7 +52,7 @@ SimulationConfig parse_command_line(int argc, char* argv[]) {
     desc.add_options()
         ("help,h", "Show help message")
         ("config,c", po::value<std::string>(&config.config_file)->required(), 
-         "Path to planet configuration file (JSON)")
+         "Path to body configuration file (JSON)")
         ("output,o", po::value<std::string>(&config.output_file)->required(), 
          "Path to output file")
         ("time,t", po::value<std::string>()->required(),
@@ -97,7 +97,7 @@ SimulationConfig parse_command_line(int argc, char* argv[]) {
     return config;
 }
 
-std::vector<Body> load_planets(const std::string& config_file) {
+std::vector<Body> load_bodies(const std::string& config_file) {
     std::ifstream f(config_file);
     if (!f.is_open()) {
         throw std::runtime_error("Could not open config file: " + config_file);
@@ -107,12 +107,12 @@ std::vector<Body> load_planets(const std::string& config_file) {
     f >> j;
     
     std::vector<Body> bodies;
-    for (const auto& planet : j["planets"]) {
-        std::string name = planet["name"];
-        double mass = planet["mass"];
+    for (const auto& body : j["bodies"]) {
+        std::string name = body["name"];
+        double mass = body["mass"];
         
-        std::vector<double> pos = planet["position"];
-        std::vector<double> vel = planet["velocity"];
+        std::vector<double> pos = body["position"];
+        std::vector<double> vel = body["velocity"];
         
         if (pos.size() != 3 || vel.size() != 3) {
             throw std::runtime_error("Position and velocity must be 3D vectors");
@@ -124,6 +124,7 @@ std::vector<Body> load_planets(const std::string& config_file) {
         Mass normalized_mass = norm_mass(mass);
         
         bodies.emplace_back(position, velocity, normalized_mass, name);
+        
         std::cout << "Loaded " << name << " with mass " << mass << " kg\n";
     }
     
@@ -135,8 +136,8 @@ int main(int argc, char* argv[]) {
         // Parse command line options
         auto config = parse_command_line(argc, argv);
         
-        // Load planets from config file
-        auto bodies = load_planets(config.config_file);
+        // Load bodies from config file
+        auto bodies = load_bodies(config.config_file);
         
         // Create and run simulator
         Simulator simulator(

src/simulator.cpp

@@ -6,16 +6,16 @@
 #include <stdexcept>
 
 Simulator::Simulator(const std::vector<Body>& bodies,
-                     Decimal step_size,
-                     int steps_per_save,
-                     const std::filesystem::path& output_file,
-                     int current_step,
-                     bool overwrite_output)
-    : bodies(bodies)
-    , step_size(step_size)
-    , steps_per_save(steps_per_save)
-    , output_file(output_file)
-    , current_step(current_step) {
+                    Decimal step_size,
+                    int steps_per_save,
+                    const std::filesystem::path& output_file,
+                    int current_step,
+                    bool overwrite_output)
+    : bodies(bodies),
+      step_size(step_size),
+      steps_per_save(steps_per_save),
+      output_file(output_file),
+      current_step(current_step) {
     
     if (std::filesystem::exists(output_file) && !overwrite_output) {
         throw std::runtime_error("File " + output_file.string() + " exists and overwrite flag not given.");
@@ -28,14 +28,14 @@ Simulator::Simulator(const std::vector<Body>& bodies,
         clear.close();
     }
     
-    // Write initial header with masses
     out.open(output_file, std::ios::app);
     if (!out) {
         throw std::runtime_error("Failed to open output file: " + output_file.string());
     }
     
     for (const auto& body : bodies) {
-        out << body.getName() << ": " << body.getMass() << "\n";
+        out << body.getName() << ": mass=" << body.getMass();
+        out << "\n";
     }
     out << "\n";
     
@@ -76,43 +76,35 @@ void Simulator::calculate_forces() {
         body.setAcceleration(Acceleration{Decimal(0), Decimal(0), Decimal(0)});
     }
     
-    for (size_t i = 0; i < bodies.size(); i++){
+    for (size_t i = 0; i < bodies.size(); i++) {
         for (size_t j = i + 1; j < bodies.size(); j++) {
             Position vec;
             for (int k = 0; k < 3; ++k) {
-                vec[k] = bodies[i].getPosition()[k] - bodies[j].getPosition()[k];
+                vec[k] = positions[i][k] - positions[j][k];
             }
             
-            // Calculate distance
             Decimal dist = std::sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
             
             // 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();
             Decimal acc_magnitude_j = force_magnitude * bodies[i].getMass();
 
             // 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);
+            bodies[i].addAcceleration(acc_i);
+            bodies[j].addAcceleration(acc_j);
         }
     }
+
 }
 
 void Simulator::move_bodies() {
@@ -125,6 +117,5 @@ void Simulator::checkpoint() {
     for (const auto& body : bodies) {
         out << body.toString() << "\n";
     }
-    out << "\n";
     out.flush();
 } 

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