Some changes - not sure exactly what

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,5 +1,5 @@
 {
-    "planets": [
+    "bodies": [
         {
             "name": "Earth",
             "mass": 5.972e24,
@@ -10,7 +10,7 @@
             "name": "Moon",
             "mass": 7.34767309e22,
             "position": [149982270700, 0, 0],
-            "velocity": [0, 30802, 0]
+            "velocity": [0, 30822, 0]
         },
         {
             "name": "Sun",

src/body.cpp

@@ -3,7 +3,6 @@
 #include <cmath>
 #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)};
@@ -15,6 +14,12 @@ void Body::addAcceleration(const Acceleration& 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);
 }
@@ -66,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 ";
+        pos_str += format_sig_figs(real_pos(X[i]), 10) + "m ";
-        vel_str += format_sig_figs(real_vel(V[i]), 3) + "m/s ";
+        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 Position& getPosition() const { return X; } 
-    const Velocity& getVelocity() const { return V; }
+    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; }
@@ -31,13 +33,13 @@ public:
     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;
+    Position X; 
-    Velocity V;
+    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"]) {
+    for (const auto& body : j["bodies"]) {
-        std::string name = planet["name"];
+        std::string name = body["name"];
-        double mass = planet["mass"];
+        double mass = body["mass"];
         
-        std::vector<double> pos = planet["position"];
+        std::vector<double> pos = body["position"];
-        std::vector<double> vel = planet["velocity"];
+        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
+        // Load bodies from config file
-        auto bodies = load_planets(config.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,
+                    Decimal step_size,
-                     int steps_per_save,
+                    int steps_per_save,
-                     const std::filesystem::path& output_file,
+                    const std::filesystem::path& output_file,
-                     int current_step,
+                    int current_step,
-                     bool overwrite_output)
+                    bool overwrite_output)
-    : bodies(bodies)
+    : bodies(bodies),
-    , step_size(step_size)
+      step_size(step_size),
-    , steps_per_save(steps_per_save)
+      steps_per_save(steps_per_save),
-    , output_file(output_file)
+      output_file(output_file),
-    , current_step(current_step) {
+      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,14 +76,13 @@ 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
@@ -91,7 +90,6 @@ void Simulator::calculate_forces() {
             // 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();
 
@@ -102,11 +100,11 @@ void Simulator::calculate_forces() {
                 acc_j[k] = vec[k] * acc_magnitude_j;   
             }
             
-            // Add this component to the two bodies
             bodies[i].addAcceleration(acc_i);
             bodies[j].addAcceleration(acc_j);
         }
     }
+
 }
 
 void Simulator::move_bodies() {
@@ -119,6 +117,5 @@ void Simulator::checkpoint() {
     for (const auto& body : bodies) {
         out << body.toString() << "\n";
     }
-    out << "\n";
     out.flush();
 } 

tests/simulator_test.cpp

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