From 460d984b3b6d8391d8e6ed30805a054055256272 Mon Sep 17 00:00:00 2001 From: Thomas Faour Date: Mon, 2 Jun 2025 23:08:35 -0400 Subject: [PATCH] Some changes - not sure exactly what --- animate.py | 36 ++++++++++++++-- config/planets.json | 4 +- src/body.cpp | 11 +++-- src/body.hpp | 14 +++--- src/main.cpp | 19 +++++---- src/simulator.cpp | 33 +++++++------- tests/simulator_test.cpp | 92 ++++++++++++++++++++-------------------- 7 files changed, 122 insertions(+), 87 deletions(-) diff --git a/animate.py b/animate.py index 2ae5a3d..0b57d29 100644 --- a/animate.py +++ b/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() \ No newline at end of file diff --git a/config/planets.json b/config/planets.json index 6e44305..8b763a7 100644 --- a/config/planets.json +++ b/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", diff --git a/src/body.cpp b/src/body.cpp index 69e7b6c..41d2069 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -3,7 +3,6 @@ #include #include - 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 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 "; - 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; } \ No newline at end of file diff --git a/src/body.hpp b/src/body.hpp index 786ffc2..69dec10 100644 --- a/src/body.hpp +++ b/src/body.hpp @@ -3,10 +3,12 @@ #include "units.hpp" #include #include +#include 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 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; } @@ -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; - Velocity V; + Position X; + Velocity V; Acceleration A; Mass m; std::string name; diff --git a/src/main.cpp b/src/main.cpp index 7ea361b..4a98530 100644 --- a/src/main.cpp +++ b/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(&config.config_file)->required(), - "Path to planet configuration file (JSON)") + "Path to body configuration file (JSON)") ("output,o", po::value(&config.output_file)->required(), "Path to output file") ("time,t", po::value()->required(), @@ -97,7 +97,7 @@ SimulationConfig parse_command_line(int argc, char* argv[]) { return config; } -std::vector load_planets(const std::string& config_file) { +std::vector 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 load_planets(const std::string& config_file) { f >> j; std::vector 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 pos = planet["position"]; - std::vector vel = planet["velocity"]; + std::vector pos = body["position"]; + std::vector 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 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( diff --git a/src/simulator.cpp b/src/simulator.cpp index 931651e..15fe6bc 100644 --- a/src/simulator.cpp +++ b/src/simulator.cpp @@ -6,16 +6,16 @@ #include Simulator::Simulator(const std::vector& 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& 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(); } \ No newline at end of file diff --git a/tests/simulator_test.cpp b/tests/simulator_test.cpp index 828bdaf..146bde2 100644 --- a/tests/simulator_test.cpp +++ b/tests/simulator_test.cpp @@ -1,54 +1,54 @@ -#include -#include "simulator.hpp" -#include -#include -#include +// #include +// #include "simulator.hpp" +// #include +// #include +// #include -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( - bodies, - Decimal(0.1), // step_size - 1, // steps_per_save - std::filesystem::path("test_output.json"), // output_file - 0, // current_step - true // overwrite_output - ); - } +// 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( +// 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 bodies; - std::unique_ptr simulator; -}; +// std::vector bodies; +// std::unique_ptr simulator; +// }; -TEST_F(SimulatorTest, Step) { - // Take a step - //TODO -} +// TEST_F(SimulatorTest, Step) { +// // Take a step +// //TODO +// } -TEST_F(SimulatorTest, TotalEnergy) { - //TODO - // Decimal initial_energy = simulator->total_energy(); +// TEST_F(SimulatorTest, TotalEnergy) { +// //TODO +// // 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); +// }