diff --git a/.gitignore b/.gitignore index bdf3094..8783d80 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,16 @@ __pycache__ /build /cmake-build-* CMakeFiles/ -*.cmake \ No newline at end of file +*.cmake +/rust/target +/rust/target/* +rust/Cargo.lock +out.out +output.json +*.png + + +# Added by cargo + +/target +Cargo.lock \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt deleted file mode 100644 index b936e89..0000000 --- a/CMakeLists.txt +++ /dev/null @@ -1,85 +0,0 @@ -cmake_minimum_required(VERSION 3.10) -project(orbital_simulator) - -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 - src/main.cpp - src/body.cpp - src/calc.cpp - 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 - src/calc.hpp - src/simulator.hpp - src/units.hpp -) - -# Create main executable -add_executable(orbital_simulator ${SOURCES} ${HEADERS}) - -# 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) - find_package(Curses REQUIRED) - target_compile_definitions(orbital_simulator PRIVATE NCURSES_ENABLED) - target_link_libraries(orbital_simulator PRIVATE ${CURSES_LIBRARIES}) -endif() \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..6e2bdc5 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,19 @@ +[package] +name = "orbital_simulator" +version = "0.1.0" +edition = "2021" +authors = ["thomas"] + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +bincode = "1.3" +clap = { version = "4.5.39", features = ["derive"] } +env_logger = "0.11.8" +glam = { version = "0.30.4", features = ["serde"] } +humantime = "2.2.0" +indicatif = "0.17.11" +log = "0.4.27" +once_cell = "1.21.3" +serde = { version = "1.0.219", features = ["derive"] } +serde_json = "1.0.140" diff --git a/README.md b/README.md deleted file mode 100644 index 14fb52c..0000000 --- a/README.md +++ /dev/null @@ -1,83 +0,0 @@ -# Orbital Simulator - -A C++ implementation of an N-body orbital simulator with high-precision calculations. - -## Dependencies - -- C++17 compatible compiler -- CMake 3.10 or higher -- Boost library (for high-precision decimal arithmetic) -- Optional: ncurses (for terminal plotting) - -## Building - -1. Create a build directory: -```bash -mkdir build -cd build -``` - -2. Configure with CMake: -```bash -cmake .. -``` - -3. Build the project: -```bash -make -``` - -To enable terminal plotting with ncurses, configure with: -```bash -cmake -DENABLE_NCURSES=ON .. -``` - -## Usage - -The simulator can be used to simulate orbital mechanics with high precision. The main components are: - -- `Body`: Represents a celestial body with position, velocity, and mass -- `Simulator`: Manages the simulation of multiple bodies -- `units.hpp`: Contains physical constants and unit conversion utilities - -Example usage: - -```cpp -#include "simulator.hpp" -#include "units.hpp" - -int main() { - // Create bodies - std::vector bodies; - - // Earth - 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"); - - // Moon - Position moon_pos{AU, Decimal(0), Decimal(0)}; - Velocity moon_vel{Decimal(0), MOON_ORBITAL_VELOCITY, Decimal(0)}; - bodies.emplace_back(moon_pos, moon_vel, MOON_MASS, "Moon"); - - // Create simulator - Simulator sim(bodies, 0.1, 100, "output.txt"); - - // Run simulation - sim.run(1000); - - return 0; -} -``` - -## Features - -- High-precision decimal arithmetic using Boost.Multiprecision -- N-body gravitational simulation -- Progress tracking and checkpointing -- Optional terminal visualization -- Configurable simulation parameters - -## License - -This project is open source and available under the MIT License. \ No newline at end of file diff --git a/animate.py b/animate.py deleted file mode 100644 index 0b57d29..0000000 --- a/animate.py +++ /dev/null @@ -1,165 +0,0 @@ -#!/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'] - - # Store position - positions[name].append((pos[0], pos[1])) - frame += 1 - times.append(frame) - - return positions, times - -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)) - - # 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)') - title = 'Orbital Simulation' - if center_body: - title += f' (centered on {center_body})' - ax.set_title(title) - 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)') - 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, args.center) - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/config/jwst.json b/config/jwst.json new file mode 100644 index 0000000..6741461 --- /dev/null +++ b/config/jwst.json @@ -0,0 +1,22 @@ +{ + "bodies": [ + { + "name": "Earth", + "mass": 5.972e24, + "position": [1.47095e11, 0, 0], + "velocity": [0, 29290, 0] + }, + { + "name": "Sun", + "mass": 1.989e30, + "position": [0, 0, 0], + "velocity": [0, 0, 0] + }, + { + "name": "JWST", + "mass": 6500, + "position": [149217067274.40607, 0, 0], + "velocity": [0, 29729.784, 0] + } + ] +} \ No newline at end of file diff --git a/src/body.cpp b/src/body.cpp deleted file mode 100644 index 41d2069..0000000 --- a/src/body.cpp +++ /dev/null @@ -1,78 +0,0 @@ -#include "body.hpp" -#include "calc.hpp" -#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)}; -} - -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 Body::save() const { - return std::make_tuple(X, V, m); -} - -Body Body::load(const std::tuple& tup) { - return Body(std::get<0>(tup), std::get<1>(tup), std::get<2>(tup)); -} - -void Body::step(Decimal step_size) { - for (int i = 0; i < 3; ++i) { - X[i] += step_size * V[i]; - V[i] += step_size * A[i]; - A[i] = Decimal(0); - } -} - -Decimal Body::E() const { - return ke() + pe(); -} - -Decimal Body::pe() const { - return -m / dist_from_o(); -} - -Decimal Body::dist_from_o() const { - Decimal sum = Decimal(0); - for (const auto& x : X) { - sum += x * x; - } - return std::sqrt(sum); -} - -Decimal Body::ke() const { - return Decimal(0.5) * m * (_speed() * _speed()); -} - -Decimal Body::_speed() const { - Decimal sum = Decimal(0); - for (const auto& v : V) { - sum += v * v; - } - return std::sqrt(sum); -} - -std::string Body::speed() const { - return format_sig_figs(real_vel(_speed()), 5); -} - -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]), 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 deleted file mode 100644 index 69dec10..0000000 --- a/src/body.hpp +++ /dev/null @@ -1,48 +0,0 @@ -#pragma once - -#include "units.hpp" -#include -#include -#include - -class Body { -public: - Body(const Position& X, const Velocity& V, const Mass& m, - const std::string& name = ""); - - // Save and load state - std::tuple save() const; - static Body load(const std::tuple& tup); - - // Physics calculations - void step(Decimal step_size); - Decimal E() const; // Total energy - Decimal pe() const; // Potential energy - Decimal ke() const; // Kinetic energy - Decimal dist_from_o() const; // Distance from origin - std::string speed() const; // Speed as formatted string - - // Getters - 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; } - - // 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; - Acceleration A; - Mass m; - std::string name; - - Decimal _speed() const; // Internal speed calculation -}; \ No newline at end of file diff --git a/src/calc.cpp b/src/calc.cpp deleted file mode 100644 index 7526c65..0000000 --- a/src/calc.cpp +++ /dev/null @@ -1,147 +0,0 @@ -#include "calc.hpp" -#include -#include -#include -#include -#include - -std::vector> calculate_distances(const std::vector& positions) { - int N = positions.size(); - std::vector> dists(N, std::vector(N, Decimal(0))); - - for (int i = 0; i < N; ++i) { - for (int j = i + 1; j < N; ++j) { - Decimal sum = Decimal(0); - for (int k = 0; k < 3; ++k) { - Decimal diff = positions[i][k] - positions[j][k]; - sum += diff * diff; - } - Decimal d = std::sqrt(sum); - dists[i][j] = d; - dists[j][i] = d; - } - } - return dists; -} - -std::string format_sig_figs(Decimal value, int sig_figs) { - if (value == 0) { - return "0"; - } - - std::stringstream ss; - ss << std::scientific << std::setprecision(sig_figs - 1) << value; - return ss.str(); -} - -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; - - std::string bar; - bar.reserve(length + 1); - bar += '['; - bar.append(filled_length, '#'); - bar.append(length - filled_length, '-'); - bar += ']'; - - auto now = std::chrono::steady_clock::now(); - auto elapsed = std::chrono::duration_cast(now - start_time).count(); - double steps_per_second = elapsed > 0 ? real_time(Decimal(iteration) * step_size) / elapsed : 0; - - // 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 -void plot_points_terminal(const std::vector& vectors, WINDOW* stdscr, - Decimal scale, int grid_width, int grid_height) { - if (vectors.empty()) { - mvwaddstr(stdscr, 0, 0, "No vectors provided."); - wrefresh(stdscr); - return; - } - - // Scale and round vectors - std::vector> scaled_vectors; - for (const auto& vec : vectors) { - scaled_vectors.emplace_back( - std::round(vec[0] / scale), - std::round(vec[1] / scale) - ); - } - - // Find bounds - int min_x = scaled_vectors[0].first; - int max_x = min_x; - int min_y = scaled_vectors[0].second; - int max_y = min_y; - - for (const auto& vec : scaled_vectors) { - min_x = std::min(min_x, vec.first); - max_x = std::max(max_x, vec.first); - min_y = std::min(min_y, vec.second); - max_y = std::max(max_y, vec.second); - } - - // Center offsets - int center_x = (grid_width / 2) - min_x; - int center_y = (grid_height / 2) - min_y; - - // Adjust coordinates - std::vector> adjusted_vectors; - for (const auto& vec : scaled_vectors) { - adjusted_vectors.emplace_back( - vec.first + center_x, - vec.second + center_y - ); - } - - // Get terminal bounds - int max_terminal_y, max_terminal_x; - getmaxyx(stdscr, max_terminal_y, max_terminal_x); - max_x = std::min(grid_width, max_terminal_x - 5); - max_y = std::min(grid_height, max_terminal_y - 5); - - // Draw grid - for (int i = grid_height; i >= 0; --i) { - std::string row = std::to_string(i - center_y) + " | "; - for (int j = 0; j <= grid_width; ++j) { - bool has_point = false; - for (const auto& vec : adjusted_vectors) { - if (vec.first == j && vec.second == i) { - has_point = true; - break; - } - } - row += has_point ? "● " : ". "; - } - mvwaddstr(stdscr, max_y - i, 0, row.substr(0, max_terminal_x - 1).c_str()); - } - - // Print X-axis labels - std::string x_labels = " "; - for (int j = 0; j <= max_x; ++j) { - x_labels += std::to_string(j - center_x) + " "; - } - mvwaddstr(stdscr, max_y + 1, 0, x_labels.substr(0, max_terminal_x - 1).c_str()); - - wrefresh(stdscr); -} -#endif \ No newline at end of file diff --git a/src/calc.hpp b/src/calc.hpp deleted file mode 100644 index 054c322..0000000 --- a/src/calc.hpp +++ /dev/null @@ -1,22 +0,0 @@ -#pragma once - -#include "units.hpp" -#include -#include -#include - -// Calculate distances between all bodies -std::vector> calculate_distances(const std::vector& positions); - -// Format a number to a specified number of significant figures -std::string format_sig_figs(Decimal value, int sig_figs); - -// Print progress bar -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 -#include -void plot_points_terminal(const std::vector& vectors, WINDOW* stdscr, - Decimal scale = 500000, int grid_width = 30, int grid_height = 30); -#endif \ No newline at end of file diff --git a/src/config.rs b/src/config.rs new file mode 100644 index 0000000..4182ec6 --- /dev/null +++ b/src/config.rs @@ -0,0 +1,19 @@ +use crate::types; +use serde::{Deserialize, Serialize}; + + +#[derive(Debug, Deserialize, Serialize)] +pub struct Body { + pub name: String, + pub mass: types::Mass, + pub position: types::Position, + pub velocity: types::Velocity, + + #[serde(default)] + pub acceleration: types::Acceleration, +} + +#[derive(Debug, Deserialize)] +pub struct ConfigFile { + pub bodies: Vec, +} \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp deleted file mode 100644 index 4a98530..0000000 --- a/src/main.cpp +++ /dev/null @@ -1,167 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include "simulator.hpp" -#include "body.hpp" -#include "units.hpp" - -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 body 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 { - po::store(po::parse_command_line(argc, argv, desc), vm); - - if (vm.count("help")) { - std::cout << desc << "\n"; - exit(0); - } - - po::notify(vm); - - // Parse simulation time - config.simulation_time = parse_time(vm["time"].as()); - - // Convert step size to Decimal and normalize - config.step_size = norm_time(Decimal(temp_step_size)); - - // 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_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); - } - - json j; - f >> j; - - std::vector bodies; - for (const auto& body : j["bodies"]) { - std::string name = body["name"]; - double mass = body["mass"]; - - 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"); - } - - // 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 bodies from config file - auto bodies = load_bodies(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; - } -} \ No newline at end of file diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..ebf6bd1 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,122 @@ +// Standard library +use std::error::Error; +use std::fs::File; +use std::io::BufReader; +use std::path::Path; +use std::time::{Duration,Instant}; + +// External crates +use clap::Parser; +use log::{info, debug}; +use indicatif::{ProgressBar, ProgressStyle}; +use humantime; + +// Internal modules +mod types; +mod config; +mod simulator; + +// Specific uses from modules +use crate::simulator::Simulation; +use crate::types::{norm_time, real_body}; + + +#[derive(Parser, Debug)] +#[command( + version, + about="Orbital mechanics simulator", + long_about = "Given initial conditions you provide to --config, \ +this program will numerically integrate and determinate their \ +paths based off Newton's law of gravity.")] +struct Args { + ///Config file for initial conditions + #[arg(short, long)] + config: String, + + /// Time to run simulation for (e.g. 10s, 5m, 2h, 100d) + #[arg(short, long, value_parser = humantime::parse_duration)] + time: Duration, + + ///Step size for simulation (seconds) + #[arg(short, long, default_value_t = 10.0)] + step_size: f64, + + ///How often to update progress bar and save (seconds) + #[arg(short = 'P', long, default_value_t = 1000)] + steps_per_save: usize, + + /// Filename to save trajectory to + #[arg(short, long)] + output_file: String, +} + +fn read_config>(path: P) -> Result> { + let file = File::open(path)?; + let reader = BufReader::new(file); + + let u = serde_json::from_reader(reader)?; + + Ok(u) +} +//fn parse_time(arg: &str) + +fn format_duration_single_unit(dur: Duration) -> String { + let secs = dur.as_secs_f64(); + + const MINUTE: f64 = 60.0; + const HOUR: f64 = 60.0 * MINUTE; + const DAY: f64 = 24.0 * HOUR; + const YEAR: f64 = 365.25 * DAY; + + if secs >= YEAR { + format!("{:.0} years", (secs / YEAR).round()) + } else if secs >= DAY { + format!("{:.0} days", (secs / DAY).round()) + } else if secs >= HOUR { + format!("{:.0} hours", (secs / HOUR).round()) + } else if secs >= MINUTE { + format!("{:.0} minutes", (secs / MINUTE).round()) + } else { + format!("{:.0} seconds", secs.round()) + } +} + +fn main() { + env_logger::init(); + let args = Args::parse(); + + info!("Loading initial parameters from {}", args.config); + + let conf = read_config(args.config).unwrap(); + for body in &conf.bodies { + info!("Loaded {} with mass {:.3e} kg", body.name, body.mass); + debug!("R_i = {:?}, V_i = {:?}", body.position, body.velocity); + } + + let n_steps = (args.time.as_secs() as f64 / args.step_size) as usize; + + + let pb = ProgressBar::new(n_steps.try_into().unwrap()); + pb.set_style( + ProgressStyle::with_template("[{elapsed_precise}] {bar:40.cyan/blue} {percent_precise}% \n\ + Time remaining: {eta} \n\ + Current simulation speed: {msg}") + .unwrap() + .progress_chars("=>-"), + ); + + let mut sim = Simulation::new( + conf.bodies, + norm_time(args.step_size), + args.steps_per_save, + args.output_file, + ); + let start = Instant::now(); + sim.run(n_steps, Some(|| { + let elapsed = start.elapsed().as_secs() as f64 + 1.0; + let speed = Duration::from_secs_f64(pb.position() as f64 * args.step_size / elapsed); + pb.set_message(format!("{} /sec", format_duration_single_unit(speed))); + pb.inc(args.steps_per_save as u64); + }) + ); +} diff --git a/src/simulator.cpp b/src/simulator.cpp deleted file mode 100644 index 15fe6bc..0000000 --- a/src/simulator.cpp +++ /dev/null @@ -1,121 +0,0 @@ -#include "simulator.hpp" -#include "calc.hpp" -#include -#include -#include -#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) { - - if (std::filesystem::exists(output_file) && !overwrite_output) { - throw std::runtime_error("File " + output_file.string() + " exists and overwrite flag not given."); - } - - if (std::filesystem::exists(output_file) && overwrite_output) { - std::cout << "Warning! Overwriting file: " << output_file.string() << std::endl; - // Clear the file if we're overwriting - std::ofstream clear(output_file, std::ios::trunc); - clear.close(); - } - - 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() << ": mass=" << body.getMass(); - out << "\n"; - } - out << "\n"; - - // Write initial state - checkpoint(); -} - -Simulator Simulator::from_checkpoint(const std::filesystem::path& output_file) { - // TODO: Implement checkpoint loading - // This would require implementing a binary format for saving/loading checkpoints - throw std::runtime_error("Checkpoint loading not implemented yet"); -} - -void Simulator::run(int steps) { - auto start_time = std::chrono::steady_clock::now(); - - for (int i = 0; i < steps; ++i) { - calculate_forces(); - move_bodies(); - if (i % steps_per_save == 0) { - checkpoint(); - print_progress_bar(i + 1, steps, start_time, 50, step_size); - } - } -} - -void Simulator::calculate_forces() { - std::vector positions; - positions.reserve(bodies.size()); - for (const auto& body : bodies) { - positions.push_back(body.getPosition()); - } - - auto dists = calculate_distances(positions); - - // Reset all accelerations to zero - for (auto& body : bodies) { - body.setAcceleration(Acceleration{Decimal(0), Decimal(0), Decimal(0)}); - } - - 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] = positions[i][k] - positions[j][k]; - } - - 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. Use r^3 to avoid normalizing vec when multiplying later - Decimal force_magnitude = 1 / (dist * dist * dist); - - 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; - acc_j[k] = vec[k] * acc_magnitude_j; - } - - bodies[i].addAcceleration(acc_i); - bodies[j].addAcceleration(acc_j); - } - } - -} - -void Simulator::move_bodies() { - for (auto& body : bodies) { - body.step(step_size); - } -} - -void Simulator::checkpoint() { - for (const auto& body : bodies) { - out << body.toString() << "\n"; - } - out.flush(); -} \ No newline at end of file diff --git a/src/simulator.hpp b/src/simulator.hpp deleted file mode 100644 index dc42243..0000000 --- a/src/simulator.hpp +++ /dev/null @@ -1,43 +0,0 @@ -#pragma once - -#include "body.hpp" -#include -#include -#include -#include -#include - -class Simulator { -public: - Simulator(const std::vector& bodies, - Decimal step_size, - int steps_per_save, - const std::filesystem::path& output_file, - int current_step = 0, - bool overwrite_output = false); - - // Create simulator from checkpoint - static Simulator from_checkpoint(const std::filesystem::path& output_file); - - // Run simulation - void run(int steps); - - // Getters - const std::vector& 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(); - void checkpoint(); - - std::vector bodies; - Decimal step_size; - int steps_per_save; - std::filesystem::path output_file; - std::ofstream out; - int current_step; -}; \ No newline at end of file diff --git a/src/simulator.rs b/src/simulator.rs new file mode 100644 index 0000000..a3f6fb5 --- /dev/null +++ b/src/simulator.rs @@ -0,0 +1,113 @@ +use std::fs::OpenOptions; +use std::io::BufWriter; + +use serde::Serialize; +use log::{debug, trace}; +use glam::DVec3; + +use crate::config::Body; +use crate::types::{norm_mass, norm_pos, norm_vel, norm_time, real_pos, real_vel, real_time, real_body}; + +pub struct Simulation { + pub bodies: Vec, + step_size: f64, + steps_per_save: usize, + save_file: String +} + +#[derive(Serialize)] +struct Snapshot<'a> { + time: f64, + bodies: &'a [Body], +} + +impl Simulation { + pub fn new( + mut bodies: Vec, + mut step_size: f64, + steps_per_save: usize, + save_file: String, + ) -> Self { + let n = bodies.len(); + for i in 0..n { + bodies[i].mass = norm_mass(bodies[i].mass); + bodies[i].position = norm_pos(bodies[i].position); + bodies[i].velocity = norm_vel(bodies[i].velocity); + } + step_size = norm_time(step_size); + Self {bodies, step_size, steps_per_save, save_file} + } + + pub fn run(&mut self, steps: usize, mut on_step: Option) + where F: FnMut() { + let file = OpenOptions::new() + .create(true) + .append(true) + .open(&self.save_file); + let mut writer = BufWriter::new(file.unwrap()); + + for i in 0..steps { + debug!("Step: {}", i); + self.reset_accelerations(); + self.caclulate_accelerations(); + self.step_bodies(); + if i % self.steps_per_save == 0 { + //save the state + let real_bodies: Vec = self.bodies.iter().map(real_body).collect(); + let snapshot = Snapshot { + time: real_time(i as f64*self.step_size), + bodies: &real_bodies, + }; + + bincode::serialize_into(&mut writer, &snapshot).expect("Couldn't write to trajectory. "); + //Do the progress bar + if let Some(f) = &mut on_step { + f(); + } + } + } + } + + fn caclulate_accelerations(&mut self) { + let r_hat_over_r3 = self.get_rhat_over_r_three(); + let n = self.bodies.len(); + for i in 0..(n-1) { + let mass_i = self.bodies[i].mass; + for j in (i+1)..n { + let mass_j = self.bodies[j].mass; + self.bodies[i].acceleration += r_hat_over_r3[i][j]*mass_j; + self.bodies[j].acceleration -= r_hat_over_r3[i][j]*mass_i; + } + } + } + + fn reset_accelerations(&mut self) { + for i in 0..self.bodies.len() { + self.bodies[i].acceleration = DVec3::ZERO; + } + } + + fn get_rhat_over_r_three (&self) -> Vec>{ + let n = self.bodies.len(); + let mut r_hat_over_r3 = vec![vec![DVec3::ZERO; n]; n]; + for i in 0..(n-1) { + for j in (i+1)..n { + let r_ij = self.bodies[j].position - self.bodies[i].position; + let dist = r_ij.length(); + r_hat_over_r3[i][j] = r_ij / dist.powf(3.0); + } + } + trace!("rhat over r3: {:?}", r_hat_over_r3); + return r_hat_over_r3; + } + + fn step_bodies(&mut self) { + for body in &mut self.bodies { + //do for each of three (x,y,z) dimensions + body.position += self.step_size*body.velocity; + body.velocity += self.step_size*body.acceleration; + trace!("{} now at {:?}", body.name, real_pos(body.position)); + trace!("{} moving at {:?}", body.name, real_vel(body.velocity)); + } + } +} \ No newline at end of file diff --git a/src/types.rs b/src/types.rs new file mode 100644 index 0000000..fa6a245 --- /dev/null +++ b/src/types.rs @@ -0,0 +1,105 @@ +#![allow(unused)] +use glam::DVec3; +use once_cell::sync::Lazy; + +use crate::config::Body; + +pub type Position = DVec3; +pub type Velocity = DVec3; +pub type Acceleration = DVec3; +pub type Mass = f64; +pub type Time = f64; + + +// Constants +const EARTH_RADIUS: f64 = 6.378e6; //meters +const EARTH_MASS: f64 = 5.972e24; //kg +const EARTH_ORBITAL_VELOCITY: f64 = 2.9780e4; //m/s +const AU: f64 = 1.49597870700e11; //meters + +const MOON_MASS: Mass = 7.34767309e22; // kg +const MOON_ORBITAL_VELOCITY: f64 = 1.022e3; //m/s relative to earth + +const SUN_MASS: f64 = 1.989e30; //kg +const SUN_RADIUS: f64 = 6.957e8; //meters + +const G: f64 = 6.67430e-11; +const R_0: f64 = EARTH_RADIUS; +const M_0: f64 = EARTH_MASS; + +static T_0: Lazy = Lazy::new(|| { + (R_0.powf(3.0) / (G * M_0)).sqrt() +}); + + +#[inline] +pub fn norm_pos(pos: Position) -> Position { + pos / R_0 +} + +#[inline] +pub fn real_pos(pos: Position) -> Position { + pos * R_0 +} + +#[inline] +pub fn norm_mass(mass: Mass) -> Mass { + mass / M_0 +} + +#[inline] +pub fn real_mass(mass: Mass) -> Mass { + mass * M_0 +} + +#[inline] +pub fn norm_time(time: Time) -> Time { + time / *T_0 +} + +#[inline] +pub fn real_time(time: Time) -> Time { + time * *T_0 +} + +#[inline] +pub fn norm_vel(vel: Velocity) -> Velocity { + vel / (R_0 / *T_0) +} + +#[inline] +pub fn real_vel(vel: Velocity) -> Velocity { + vel * (R_0 / *T_0) +} + +#[inline] +pub fn norm_acc(acc: Acceleration) -> Acceleration { + acc / (R_0 / (*T_0 * *T_0)) +} + +#[inline] +pub fn real_acc(acc: Acceleration) -> Acceleration { + acc * (R_0 / (*T_0 * *T_0)) +} + +#[inline] +pub fn norm_body(body: Body) -> Body { + Body { + name: body.name, + mass: norm_mass(body.mass), + position: norm_pos(body.position), + velocity: norm_vel(body.velocity), + acceleration: DVec3::ZERO, + } +} + +#[inline] +pub fn real_body(body: &Body) -> Body { + Body { + name: body.name.clone(), + mass: real_mass(body.mass), + position: real_pos(body.position), + velocity: real_vel(body.velocity), + acceleration: DVec3::ZERO, + } +} \ No newline at end of file diff --git a/src/units.hpp b/src/units.hpp deleted file mode 100644 index b9905d5..0000000 --- a/src/units.hpp +++ /dev/null @@ -1,46 +0,0 @@ -#pragma once - -#include -#include -#include -//#include - -using Decimal = long double; //boost::multiprecision::cpp_dec_float_50; - -// Type aliases for clarity -using Position = std::array; -using Velocity = std::array; -using Acceleration = std::array; -using Mass = Decimal; - -// Constants -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 = 734767309e14;//Decimal("734767309e14"); -const Decimal MOON_ORBITAL_VELOCITY = 1022;//Decimal("1022"); // m/s relative to earth - -const Decimal SUN_MASS = 1989e27;//Decimal("1989e27"); // kg -const Decimal SUN_RADIUS = 6957e5;//Decimal("6957e5"); // meters - -const Decimal PI = 3.14159265358979323846264338327950288419716939937510; - -// Normalizing constants -const Decimal G = 6.67430e-11; -const Decimal r_0 = EARTH_RADIUS; -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; } -inline Decimal real_pos(Decimal pos) { return pos * r_0; } -inline Decimal norm_mass(Decimal mass) { return mass / m_0; } -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 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)); } \ No newline at end of file diff --git a/test.py b/test.py deleted file mode 100644 index e76470a..0000000 --- a/test.py +++ /dev/null @@ -1,44 +0,0 @@ -from orbiter.orbits.body import Body -from orbiter.orbits.simulator import Simulator -from orbiter.units import * - -from pathlib import Path -from decimal import Decimal, getcontext - -getcontext().prec = 50 - -#set up the earth -earth = Body( - Position([0,0,0]), - Velocity([0,0,0]), - Mass(norm_mass(EARTH_MASS)), - "Earth" -) - -r = EARTH_RADIUS+100_000 -#Lets try a body just outside earth accelerating in. Should be 9.8m/s2 -person = Body( - Position([norm_pos(r),0,0]), #10_000m in the sky, airliner height! - Velocity([0,(Decimal(0.5)/norm_pos(r)).sqrt(),0]), #orbital velocity - Mass(norm_mass(80)), #avg person - "Person" -) - -T = 2*pi_approx*norm_pos(r)/person.V[1] - -time_to_run = T*3 #norm_time(2000) -STEP_SIZE = Decimal(6e-4) -n_steps = int(time_to_run/STEP_SIZE) - -def main(): - print("Before: ") - print(str(person)) - print(str(earth)) - import cProfile - s = Simulator([earth,person], STEP_SIZE, 100, Path("hello_world")) - cProfile.run(s.run(n_steps)) - print("\nAfter:") - print(str(person)) - print(str(earth)) - -main() \ No newline at end of file diff --git a/tests/body_test.cpp b/tests/body_test.cpp deleted file mode 100644 index 09a1109..0000000 --- a/tests/body_test.cpp +++ /dev/null @@ -1,111 +0,0 @@ -#include -#include "body.hpp" -#include -#include - -class BodyTest : public testing::Test { -protected: - BodyTest() { - test_body = std::make_unique( - Position{Decimal(0), Decimal(0), Decimal(0)}, - Velocity{Decimal(0), Decimal(0), Decimal(0)}, - Mass(1.0), - "test_body" - ); - } - - std::unique_ptr 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::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::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); -} \ No newline at end of file diff --git a/tests/calc_test.cpp b/tests/calc_test.cpp deleted file mode 100644 index fe889b5..0000000 --- a/tests/calc_test.cpp +++ /dev/null @@ -1,10 +0,0 @@ -#include -#include "calc.hpp" -#include - -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"); -} diff --git a/tests/simulator_test.cpp b/tests/simulator_test.cpp deleted file mode 100644 index 146bde2..0000000 --- a/tests/simulator_test.cpp +++ /dev/null @@ -1,54 +0,0 @@ -// #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 -// ); -// } - -// std::vector bodies; -// std::unique_ptr 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); -// } - diff --git a/tracjtory.out b/tracjtory.out new file mode 100644 index 0000000..e1a0233 Binary files /dev/null and b/tracjtory.out differ diff --git a/tracjtory.out2 b/tracjtory.out2 new file mode 100644 index 0000000..0e1ce6a Binary files /dev/null and b/tracjtory.out2 differ