diff --git a/CMakeLists.txt b/CMakeLists.txt index e18eabe..8dc063e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,8 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) # Find required packages -find_package(Boost REQUIRED COMPONENTS system) +find_package(Boost REQUIRED COMPONENTS program_options) +find_package(nlohmann_json REQUIRED) # Add source files set(SOURCES @@ -27,15 +28,16 @@ set(HEADERS add_executable(orbital_simulator ${SOURCES} ${HEADERS}) # Link libraries -target_link_libraries(orbital_simulator PRIVATE - Boost::system - ${Boost_LIBRARIES} +target_link_libraries(orbital_simulator + PRIVATE + Boost::program_options + nlohmann_json::nlohmann_json ) # Include directories -target_include_directories(orbital_simulator PRIVATE - ${CMAKE_SOURCE_DIR} - ${Boost_INCLUDE_DIRS} +target_include_directories(orbital_simulator + PRIVATE + ${CMAKE_SOURCE_DIR}/src ) # Optional: Enable ncurses for terminal plotting diff --git a/animate.py b/animate.py new file mode 100644 index 0000000..2ae5a3d --- /dev/null +++ b/animate.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python3 + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation +from collections import defaultdict +import argparse +import re + +def parse_line(line): + """Parse a line from the simulation output.""" + # Skip comments and empty lines + if line.startswith('#') or not line.strip(): + return None + + # Parse the line format: body_name: X = x1m x2m x3m, V = v1m/s v2m/s v3m/s + try: + # Split on colon to separate name and data + name_part, data_part = line.strip().split(':') + name = name_part.strip() + + # Split data part into position and velocity + pos_part, vel_part = data_part.split(',') + + # Extract position values + pos_str = pos_part.split('=')[1].strip() + pos_values = [float(x.replace('m', '').strip()) for x in pos_str.split() if x.strip()] + + if len(pos_values) != 3: + return None + + return { + 'name': name, + 'position': tuple(pos_values) + } + except (ValueError, IndexError): + return None + +def read_output_file(filename): + """Read the simulation output file and organize data by body and time.""" + positions = defaultdict(list) + times = [] # We'll use frame numbers as times since actual time isn't in output + + with open(filename, 'r') as f: + frame = 0 + for line in f: + data = parse_line(line) + if data is None: + continue + + name = data['name'] + pos = data['position'] + positions[name].append((pos[0], pos[1])) + frame += 1 + times.append(frame) + + return positions, times + +def create_animation(positions, times, output_file=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 + + # Set up the figure and axis + fig, ax = plt.subplots(figsize=(10, 10)) + + # Set equal aspect ratio + ax.set_aspect('equal') + + # Create scatter plots for each body + scatters = {} + for name in positions.keys(): + scatters[name], = ax.plot([], [], 'o-', label=name, alpha=0.7) + + # Set up the plot + ax.set_xlabel('X (m)') + ax.set_ylabel('Y (m)') + ax.set_title('Orbital Simulation') + ax.legend() + + # Find the bounds of the plot + all_x = [] + all_y = [] + for pos_list in positions.values(): + if pos_list: # Only process if we have positions + x, y = zip(*pos_list) + all_x.extend(x) + all_y.extend(y) + + if not all_x or not all_y: + print("Error: No valid position data found") + return + + max_range = max(max(abs(min(all_x)), abs(max(all_x))), + max(abs(min(all_y)), abs(max(all_y)))) + ax.set_xlim(-max_range, max_range) + ax.set_ylim(-max_range, max_range) + + def init(): + """Initialize the animation.""" + for scatter in scatters.values(): + scatter.set_data([], []) + return list(scatters.values()) + + def update(frame): + """Update the animation for each frame.""" + for name, scatter in scatters.items(): + if positions[name]: # Only update if we have positions + x, y = zip(*positions[name][:frame+1]) + scatter.set_data(x, y) + return list(scatters.values()) + + # Create the animation + anim = FuncAnimation(fig, update, frames=len(times), + init_func=init, blit=True, + interval=50) # 50ms between frames + + if output_file: + anim.save(output_file, writer='ffmpeg', fps=30) + else: + plt.show() + +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)') + args = parser.parse_args() + + positions, times = read_output_file(args.input_file) + create_animation(positions, times, args.output) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/config/planets.json b/config/planets.json new file mode 100644 index 0000000..6e44305 --- /dev/null +++ b/config/planets.json @@ -0,0 +1,22 @@ +{ + "planets": [ + { + "name": "Earth", + "mass": 5.972e24, + "position": [149597870700, 0, 0], + "velocity": [0, 29780, 0] + }, + { + "name": "Moon", + "mass": 7.34767309e22, + "position": [149982270700, 0, 0], + "velocity": [0, 30802, 0] + }, + { + "name": "Sun", + "mass": 1.989e30, + "position": [0, 0, 0], + "velocity": [0, 0, 0] + } + ] +} \ No newline at end of file diff --git a/config/planets_fake.json b/config/planets_fake.json new file mode 100644 index 0000000..d16adeb --- /dev/null +++ b/config/planets_fake.json @@ -0,0 +1,16 @@ +{ + "planets": [ + { + "name": "Earth", + "mass": 7.342e24, + "position": [0, 0, 0], + "velocity": [0, -1022, 0] + }, + { + "name": "Moon", + "mass": 7.34767309e24, + "position": [384400000, 0, 0], + "velocity": [0, 1022, 0] + } + ] +} diff --git a/orbiter/orbits/__pycache__/body.cpython-312.pyc b/orbiter/orbits/__pycache__/body.cpython-312.pyc index 0ab8563..9991296 100644 Binary files a/orbiter/orbits/__pycache__/body.cpython-312.pyc and b/orbiter/orbits/__pycache__/body.cpython-312.pyc differ diff --git a/orbiter/orbits/__pycache__/simulator.cpython-312.pyc b/orbiter/orbits/__pycache__/simulator.cpython-312.pyc index 5749f93..1fcb39f 100644 Binary files a/orbiter/orbits/__pycache__/simulator.cpython-312.pyc and b/orbiter/orbits/__pycache__/simulator.cpython-312.pyc differ diff --git a/src/body.cpp b/src/body.cpp index 5781641..9d5dc9a 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -37,7 +37,7 @@ Decimal Body::dist_from_o() const { for (const auto& x : X) { sum += x * x; } - return boost::multiprecision::sqrt(sum); + return std::sqrt(sum); } Decimal Body::ke() const { @@ -49,7 +49,7 @@ Decimal Body::_speed() const { for (const auto& v : V) { sum += v * v; } - return boost::multiprecision::sqrt(sum); + return std::sqrt(sum); } std::string Body::speed() const { diff --git a/src/calc.cpp b/src/calc.cpp index 37b305b..7526c65 100644 --- a/src/calc.cpp +++ b/src/calc.cpp @@ -16,7 +16,7 @@ std::vector> calculate_distances(const std::vector start_time, int length) { +void print_progress_bar(int iteration, int total, std::chrono::time_point start_time, int length, Decimal step_size) { float percent = (float)iteration / total * 100; int filled_length = length * iteration / total; @@ -47,10 +47,26 @@ void print_progress_bar(int iteration, int total, std::chrono::time_point(now - start_time).count(); - int steps_per_second = elapsed > 0 ? iteration / elapsed : 0; + double steps_per_second = elapsed > 0 ? real_time(Decimal(iteration) * step_size) / elapsed : 0; - std::cout << "\r" << bar << " " << std::fixed << std::setprecision(2) - << percent << "% " << steps_per_second << " steps/s" << std::flush; + // Determine appropriate time unit + std::string time_unit; + if (steps_per_second >= 3600) { + time_unit = "hour/s"; + steps_per_second /= 3600; + } else if (steps_per_second >= 60) { + time_unit = "min/s"; + steps_per_second /= 60; + } else { + time_unit = "s/s"; + } + + // Clear the current line and move cursor to start + std::cout << "\r\033[K"; + + // Print the progress bar + std::cout << bar << " " << std::fixed << std::setprecision(2) + << percent << "% " << std::setprecision(1) << steps_per_second << " " << time_unit << std::flush; } #ifdef NCURSES_ENABLED diff --git a/src/calc.hpp b/src/calc.hpp index ddcad40..054c322 100644 --- a/src/calc.hpp +++ b/src/calc.hpp @@ -12,7 +12,7 @@ std::vector> calculate_distances(const std::vector start_time, int length = 50); +void print_progress_bar(int iteration, int total, std::chrono::time_point start_time, int length, Decimal step_size); // Terminal plotting functions (if needed) #ifdef NCURSES_ENABLED diff --git a/src/main.cpp b/src/main.cpp index fb1edcc..7ea361b 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,41 +1,164 @@ -#include "simulator.hpp" -#include "units.hpp" #include +#include +#include +#include +#include +#include +#include "simulator.hpp" +#include "body.hpp" +#include "units.hpp" -int main() { +using json = nlohmann::json; +namespace po = boost::program_options; + +struct SimulationConfig { + std::string config_file; + std::string output_file; + int steps; + int steps_per_save; + Decimal step_size; // Changed to Decimal for normalized time + bool overwrite_output; + double simulation_time; // in seconds (real time) +}; + +// Convert time string to seconds +double parse_time(const std::string& time_str) { + std::string value_str = time_str; + std::string unit; + + // Extract the unit (last character) + if (!time_str.empty()) { + unit = time_str.back(); + value_str = time_str.substr(0, time_str.length() - 1); + } + + double value = std::stod(value_str); + + // Convert to seconds based on unit + switch (unit[0]) { + case 's': return value; // seconds + case 'm': return value * 60; // minutes + case 'h': return value * 3600; // hours + case 'd': return value * 86400; // days + default: throw std::runtime_error("Invalid time unit. Use s/m/h/d for seconds/minutes/hours/days"); + } +} + +SimulationConfig parse_command_line(int argc, char* argv[]) { + SimulationConfig config; + double temp_step_size; // Temporary variable for parsing + + po::options_description desc("Orbital Simulator Options"); + desc.add_options() + ("help,h", "Show help message") + ("config,c", po::value(&config.config_file)->required(), + "Path to planet configuration file (JSON)") + ("output,o", po::value(&config.output_file)->required(), + "Path to output file") + ("time,t", po::value()->required(), + "Simulation time with unit (e.g., 1h for 1 hour, 30m for 30 minutes, 2d for 2 days)") + ("step-size,s", po::value(&temp_step_size)->default_value(1.0), + "Simulation step size in seconds") + ("steps-per-save,p", po::value(&config.steps_per_save)->default_value(100), + "Number of steps between saves") + ("overwrite,w", po::bool_switch(&config.overwrite_output), + "Overwrite output file if it exists"); + + po::variables_map vm; try { - // Create bodies - std::vector bodies; + po::store(po::parse_command_line(argc, argv, desc), vm); - // Earth at origin - Position earth_pos{Decimal(0), Decimal(0), Decimal(0)}; - Velocity earth_vel{Decimal(0), Decimal(0), Decimal(0)}; - bodies.emplace_back(earth_pos, earth_vel, EARTH_MASS, "Earth"); + if (vm.count("help")) { + std::cout << desc << "\n"; + exit(0); + } - // Moon at average distance (384,400 km) with orbital velocity - // Using average orbital velocity of 1.022 km/s - Position moon_pos{Decimal("384400000"), Decimal(0), Decimal(0)}; // 384,400 km in meters - Velocity moon_vel{Decimal(0), Decimal("1022"), Decimal(0)}; // 1.022 km/s in m/s - bodies.emplace_back(moon_pos, moon_vel, MOON_MASS, "Moon"); + po::notify(vm); - // Create simulator - // Step size: 100 seconds (small enough for accuracy) - // Steps per save: 1000 (save every ~27.8 hours) - // Total steps: 27.3 days * 24 hours * 3600 seconds / 100 seconds per step - Decimal step_size = Decimal("1"); // 100 seconds per step - int steps_per_save = 1000; // Save every 1000 steps - int total_steps = 27.3 * 24 * 3600 / 1; // Steps for one orbit + // Parse simulation time + config.simulation_time = parse_time(vm["time"].as()); - Simulator sim(bodies, step_size, steps_per_save, std::filesystem::path("output.txt")); + // Convert step size to Decimal and normalize + config.step_size = norm_time(Decimal(temp_step_size)); - // Run simulation - std::cout << "Starting simulation of one lunar orbit...\n"; - std::cout << "Step size: " << step_size << " seconds\n"; - std::cout << "Total steps: " << total_steps << " (approximately 27.3 days)\n"; - sim.run(total_steps); - std::cout << "\nSimulation complete!\n"; + // Calculate number of steps based on normalized time and step size + config.steps = static_cast(norm_time(config.simulation_time) / config.step_size); + } catch (const po::error& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << desc << "\n"; + exit(1); + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << desc << "\n"; + exit(1); + } + + return config; +} + +std::vector load_planets(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); + } + + json j; + f >> j; + + std::vector bodies; + for (const auto& planet : j["planets"]) { + std::string name = planet["name"]; + double mass = planet["mass"]; + + std::vector pos = planet["position"]; + std::vector vel = planet["velocity"]; + + if (pos.size() != 3 || vel.size() != 3) { + throw std::runtime_error("Position and velocity must be 3D vectors"); + } + + // Normalize units before creating the body + Position position{norm_pos(pos[0]), norm_pos(pos[1]), norm_pos(pos[2])}; + Velocity velocity{norm_vel(vel[0]), norm_vel(vel[1]), norm_vel(vel[2])}; + Mass normalized_mass = norm_mass(mass); + + bodies.emplace_back(position, velocity, normalized_mass, name); + std::cout << "Loaded " << name << " with mass " << mass << " kg\n"; + } + + return bodies; +} + +int main(int argc, char* argv[]) { + try { + // Parse command line options + auto config = parse_command_line(argc, argv); + + // Load planets from config file + auto bodies = load_planets(config.config_file); + + // Create and run simulator + Simulator simulator( + bodies, + config.step_size, + config.steps_per_save, + config.output_file, + 0, // current_step + config.overwrite_output + ); + + std::cout << "Starting simulation with " << bodies.size() << " bodies\n"; + std::cout << "Step size: " << real_time(config.step_size) << " seconds\n"; + std::cout << "Simulation time: " << config.simulation_time << " seconds\n"; + std::cout << "Total steps: " << config.steps << "\n"; + std::cout << "Steps per save: " << config.steps_per_save << "\n"; + + simulator.run(config.steps); + + std::cout << "Simulation completed successfully\n"; return 0; + } catch (const std::exception& e) { std::cerr << "Error: " << e.what() << std::endl; return 1; diff --git a/src/simulator.cpp b/src/simulator.cpp index 67010f3..42ccb4d 100644 --- a/src/simulator.cpp +++ b/src/simulator.cpp @@ -29,7 +29,7 @@ Simulator::Simulator(const std::vector& bodies, } // Write initial header with masses - out = std::ofstream(output_file, std::ios::app); + out.open(output_file, std::ios::app); if (!out) { throw std::runtime_error("Failed to open output file: " + output_file.string()); } @@ -55,11 +55,9 @@ void Simulator::run(int steps) { for (int i = 0; i < steps; ++i) { calculate_forces(); move_bodies(); - current_step++; - - if (current_step % steps_per_save == 0) { - print_progress_bar(i, steps, start_time); + if (i % steps_per_save == 0) { checkpoint(); + print_progress_bar(i + 1, steps, start_time, 50, step_size); } } } @@ -78,7 +76,7 @@ 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) { @@ -86,22 +84,22 @@ void Simulator::calculate_forces() { } // Calculate distance - Decimal dist = boost::multiprecision::sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); + 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 - Decimal force_magnitude = G * bodies[i].getMass() * bodies[j].getMass() / (dist * dist); + // 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); // Calculate acceleration for both bodies - // a = F/m - Decimal acc_magnitude_i = force_magnitude / bodies[i].getMass(); - Decimal acc_magnitude_j = force_magnitude / bodies[j].getMass(); - + 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; // Negative because force is attractive - acc_j[k] = vec[k] * acc_magnitude_j / dist; // Positive because force is attractive + acc_i[k] = -vec[k] * acc_magnitude_i / dist; + acc_j[k] = vec[k] * acc_magnitude_j / dist; } // Add to current accelerations @@ -124,26 +122,9 @@ void Simulator::move_bodies() { } void Simulator::checkpoint() { - // Open file in append mode - std::ofstream out(output_file, std::ios::app); - if (!out) { - throw std::runtime_error("Failed to open output file: " + output_file.string()); - } - - // Write timestamp and step number - out << "Time: " << (current_step * step_size) << "s\n"; - - // Write positions and velocities for each body for (const auto& body : bodies) { - out << body.getName() << " : "; - out << " "; - for (int i = 0; i < 3; ++i) { - out << real_pos(body.getPosition()[i]) << "m "; - } - out << " : "; - for (int i = 0; i < 3; ++i) { - out << real_vel(body.getVelocity()[i]) << "m/s "; - } - out << "\n"; + out << body.toString() << "\n"; } + out << "\n"; + out.flush(); } \ No newline at end of file diff --git a/src/simulator.hpp b/src/simulator.hpp index 01ee2d0..58a8aeb 100644 --- a/src/simulator.hpp +++ b/src/simulator.hpp @@ -5,6 +5,7 @@ #include #include #include +#include class Simulator { public: diff --git a/src/units.hpp b/src/units.hpp index bb40a7f..2097234 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -14,24 +14,24 @@ using Acceleration = std::array; using Mass = Decimal; // Constants -const Decimal EARTH_MASS = Decimal("5972e21"); // kg -const Decimal EARTH_RADIUS = Decimal("6378e3"); // meters -const Decimal EARTH_ORBITAL_VELOCITY = Decimal("29780"); // m/s -const Decimal AU = Decimal("149597870700"); // meters +const Decimal EARTH_MASS = 5972e21;//Decimal("5972e21"); // kg +const Decimal EARTH_RADIUS = 6378e3;//Decimal("6378e3"); // meters +const Decimal EARTH_ORBITAL_VELOCITY = 29780;//Decimal("29780"); // m/s +const Decimal AU = 149597870700;//Decimal("149597870700"); // meters -const Decimal MOON_MASS = Decimal("734767309e14"); -const Decimal MOON_ORBITAL_VELOCITY = Decimal("1022"); // m/s relative to earth +const Decimal MOON_MASS = 734767309e14;//Decimal("734767309e14"); +const Decimal MOON_ORBITAL_VELOCITY = 1022;//Decimal("1022"); // m/s relative to earth -const Decimal SUN_MASS = Decimal("1989e27"); // kg -const Decimal SUN_RADIUS = Decimal("6957e5"); // meters +const Decimal SUN_MASS = 1989e27;//Decimal("1989e27"); // kg +const Decimal SUN_RADIUS = 6957e5;//Decimal("6957e5"); // meters -const Decimal PI = Decimal("3.14159265358979323846264338327950288419716939937510"); +const Decimal PI = 3.14159265358979323846264338327950288419716939937510; // Normalizing constants -const Decimal G = Decimal("6.67430e-11"); +const Decimal G = 6.67430e-11; const Decimal r_0 = EARTH_RADIUS; -const Decimal m_0 = Decimal("5.972e24"); -const Decimal t_0 = boost::multiprecision::sqrt((r_0 * r_0 * r_0) / (G * m_0)); +const Decimal m_0 = 5.972e24; +const Decimal t_0 = std::sqrt((r_0 * r_0 * r_0) / (G * m_0)); // Utility functions inline Decimal norm_pos(Decimal pos) { return pos / r_0; }