24 changed files with 1225 additions and 413 deletions
--- a/.gitignore
+++ b/.gitignore
@ -4,16 +4,4 @@ __pycache__
 /build
 /cmake-build-*
 CMakeFiles/
-*.cmake
+*.cmake
 /rust/target
 /rust/target/*
 rust/Cargo.lock
 out.out
 output.json
 *.png
 # Added by cargo
 /target
 Cargo.lock
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -0,0 +1,85 @@
 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() 
--- a/Cargo.toml
+++ b/Cargo.toml
@ -1,19 +0,0 @@
 [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"
--- a/README.md
+++ b/README.md
@ -0,0 +1,83 @@
 # 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<Body> 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.
--- a/animate.py
+++ b/animate.py
@ -0,0 +1,165 @@
 #!/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() 
--- a/config/jwst.json
+++ b/config/jwst.json
@ -1,22 +0,0 @@
 {
    "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]
        }
    ]
 } 
--- a/src/body.cpp
+++ b/src/body.cpp
@ -0,0 +1,78 @@
 #include "body.hpp"
 #include "calc.hpp"
 #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)};
 }
 void Body::addAcceleration(const Acceleration& new_A) {
    for (const auto& [new_a, cur_a]: std::views::zip(new_A, this->A)){
        cur_a += new_a;
    }
 }
 void Body::subAcceleration(const Acceleration& new_A) {
    for (const auto& [new_a, cur_a]: std::views::zip(new_A, this->A)){
        cur_a -= new_a;
    }
 }
 std::tuple<Position, Velocity, Mass> Body::save() const {
    return std::make_tuple(X, V, m);
 }
 Body Body::load(const std::tuple<Position, Velocity, Mass>& 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;
 } 
--- a/src/body.hpp
+++ b/src/body.hpp
@ -0,0 +1,48 @@
 #pragma once
 #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 = "");
    // Save and load state
    std::tuple<Position, Velocity, Mass> save() const;
    static Body load(const std::tuple<Position, Velocity, Mass>& 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
 }; 
--- a/src/calc.cpp
+++ b/src/calc.cpp
@ -0,0 +1,147 @@
 #include "calc.hpp"
 #include <cmath>
 #include <iostream>
 #include <iomanip>
 #include <sstream>
 #include <boost/multiprecision/cpp_dec_float.hpp>
 std::vector<std::vector<Decimal>> calculate_distances(const std::vector<Position>& positions) {
    int N = positions.size();
    std::vector<std::vector<Decimal>> dists(N, std::vector<Decimal>(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<std::chrono::steady_clock> 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<std::chrono::seconds>(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<Position>& 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<std::pair<int, int>> 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<std::pair<int, int>> 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 
--- a/src/calc.hpp
+++ b/src/calc.hpp
@ -0,0 +1,22 @@
 #pragma once
 #include "units.hpp"
 #include <vector>
 #include <string>
 #include <chrono>
 // Calculate distances between all bodies
 std::vector<std::vector<Decimal>> calculate_distances(const std::vector<Position>& 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<std::chrono::steady_clock> start_time, int length, Decimal step_size);
 // Terminal plotting functions (if needed)
 #ifdef NCURSES_ENABLED
 #include <ncurses.h>
 void plot_points_terminal(const std::vector<Position>& vectors, WINDOW* stdscr, 
                         Decimal scale = 500000, int grid_width = 30, int grid_height = 30);
 #endif 
--- a/src/config.rs
+++ b/src/config.rs
@ -1,19 +0,0 @@
 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<Body>,
 }
--- a/src/main.cpp
+++ b/src/main.cpp
@ -0,0 +1,167 @@
 #include <iostream>
 #include <fstream>
 #include <string>
 #include <vector>
 #include <boost/program_options.hpp>
 #include <nlohmann/json.hpp>
 #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<std::string>(&config.config_file)->required(), 
         "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(),
         "Simulation time with unit (e.g., 1h for 1 hour, 30m for 30 minutes, 2d for 2 days)")
        ("step-size,s", po::value<double>(&temp_step_size)->default_value(1.0), 
         "Simulation step size in seconds")
        ("steps-per-save,p", po::value<int>(&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<std::string>());
        // 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<int>(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<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);
    }
    json j;
    f >> j;
    std::vector<Body> bodies;
    for (const auto& body : j["bodies"]) {
        std::string name = body["name"];
        double mass = body["mass"];
        std::vector<double> pos = body["position"];
        std::vector<double> vel = body["velocity"];
        if (pos.size() != 3 || vel.size() != 3) {
            throw std::runtime_error("Position and velocity must be 3D vectors");
        }
        // 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;
    }
 } 
--- a/src/main.rs
+++ b/src/main.rs
@ -1,122 +0,0 @@
 // 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<P: AsRef<Path>>(path: P) -> Result<config::ConfigFile, Box<dyn Error>> {
    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);
        })
    );
 }
--- a/src/simulator.cpp
+++ b/src/simulator.cpp
@ -0,0 +1,121 @@
 #include "simulator.hpp"
 #include "calc.hpp"
 #include <fstream>
 #include <chrono>
 #include <iostream>
 #include <stdexcept>
 Simulator::Simulator(const std::vector<Body>& bodies,
                    Decimal step_size,
                    int steps_per_save,
                    const std::filesystem::path& output_file,
                    int current_step,
                    bool overwrite_output)
    : bodies(bodies),
      step_size(step_size),
      steps_per_save(steps_per_save),
      output_file(output_file),
      current_step(current_step) {
    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<Position> 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();
 } 
--- a/src/simulator.hpp
+++ b/src/simulator.hpp
@ -0,0 +1,43 @@
 #pragma once
 #include "body.hpp"
 #include <vector>
 #include <string>
 #include <filesystem>
 #include <memory>
 #include <fstream>
 class Simulator {
 public:
    Simulator(const std::vector<Body>& 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<Body>& getBodies() const { return bodies; }
    Decimal getStepSize() const { return step_size; }
    int getStepsPerSave() const { return steps_per_save; }
    const std::filesystem::path& getOutputFile() const { return output_file; }
    int getCurrentStep() const { return current_step; }
 private:
    void calculate_forces();
    void move_bodies();
    void checkpoint();
    std::vector<Body> bodies;
    Decimal step_size;
    int steps_per_save;
    std::filesystem::path output_file;
    std::ofstream out;
    int current_step;
 }; 
--- a/src/simulator.rs
+++ b/src/simulator.rs
@ -1,113 +0,0 @@
 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<Body>,
    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<Body>, 
        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<F>(&mut self, steps: usize, mut on_step: Option<F>)
    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<Body> = 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<Vec<DVec3>>{
        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));
        }
    }
 }
--- a/src/types.rs
+++ b/src/types.rs
@ -1,105 +0,0 @@
 #![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<f64> = 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,
    }
 }
--- a/src/units.hpp
+++ b/src/units.hpp
@ -0,0 +1,46 @@
 #pragma once
 #include <array>
 #include <string>
 #include <cmath>
 //#include <boost/multiprecision/cpp_dec_float.hpp>
 using Decimal = long double; //boost::multiprecision::cpp_dec_float_50;
 // Type aliases for clarity
 using Position = std::array<Decimal, 3>;
 using Velocity = std::array<Decimal, 3>;
 using Acceleration = std::array<Decimal, 3>;
 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)); } 
--- a/test.py
+++ b/test.py
@ -0,0 +1,44 @@
 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()
--- a/tests/body_test.cpp
+++ b/tests/body_test.cpp
@ -0,0 +1,111 @@
 #include <gtest/gtest.h>
 #include "body.hpp"
 #include <cmath>
 #include <memory>
 class BodyTest : public testing::Test {
 protected:
    BodyTest() {
        test_body = std::make_unique<Body>(
            Position{Decimal(0), Decimal(0), Decimal(0)},
            Velocity{Decimal(0), Decimal(0), Decimal(0)},
            Mass(1.0),
            "test_body"
        );
    }
    std::unique_ptr<Body> test_body;
 };
 TEST_F(BodyTest, AddAcceleration) {
    Acceleration new_A{Decimal(1.0), Decimal(2.0), Decimal(3.0)};
    test_body->addAcceleration(new_A);
    const auto& A = test_body->getAcceleration();
    EXPECT_DOUBLE_EQ(A[0], 1.0);
    EXPECT_DOUBLE_EQ(A[1], 2.0);
    EXPECT_DOUBLE_EQ(A[2], 3.0);
    // Add another acceleration
    Acceleration another_A{Decimal(2.0), Decimal(3.0), Decimal(4.0)};
    test_body->addAcceleration(another_A);
    EXPECT_DOUBLE_EQ(A[0], 3.0);
    EXPECT_DOUBLE_EQ(A[1], 5.0);
    EXPECT_DOUBLE_EQ(A[2], 7.0);
 }
 TEST_F(BodyTest, Step) {
    // Set initial velocity and acceleration
    test_body->setAcceleration(Acceleration{Decimal(1.0), Decimal(2.0), Decimal(3.0)});
    // Take a step
    test_body->step(1.0);
    // Check new position
    const auto& X = test_body->getPosition();
    EXPECT_DOUBLE_EQ(X[0], 0.0);
    EXPECT_DOUBLE_EQ(X[1], 0.0);
    EXPECT_DOUBLE_EQ(X[2], 0.0);
    // Check new velocity
    const auto& V = test_body->getVelocity();
    EXPECT_DOUBLE_EQ(V[0], 1.0);
    EXPECT_DOUBLE_EQ(V[1], 2.0);
    EXPECT_DOUBLE_EQ(V[2], 3.0);
    // Check acceleration is reset
    const auto& A = test_body->getAcceleration();
    EXPECT_DOUBLE_EQ(A[0], 0.0);
    EXPECT_DOUBLE_EQ(A[1], 0.0);
    EXPECT_DOUBLE_EQ(A[2], 0.0);
 }
 TEST_F(BodyTest, Energy) {
    // Set position and velocity through save/load
    auto state = std::make_tuple(
        Position{Decimal(1.0), Decimal(0.0), Decimal(0.0)},
        Velocity{Decimal(1.0), Decimal(0.0), Decimal(0.0)},
        Mass(1.0)
    );
    test_body = std::make_unique<Body>(Body::load(state));
    // Calculate expected values
    Decimal expected_ke = Decimal(0.5); // 0.5 * m * v^2
    Decimal expected_pe = Decimal(-1.0); // -m/r
    Decimal expected_total = expected_ke + expected_pe;
    EXPECT_DOUBLE_EQ(test_body->ke(), expected_ke);
    EXPECT_DOUBLE_EQ(test_body->pe(), expected_pe);
    EXPECT_DOUBLE_EQ(test_body->E(), expected_total);
 }
 TEST_F(BodyTest, SaveAndLoad) {
    // Set some values through save/load
    auto state = std::make_tuple(
        Position{Decimal(1.0), Decimal(2.0), Decimal(3.0)},
        Velocity{Decimal(4.0), Decimal(5.0), Decimal(6.0)},
        Mass(7.0)
    );
    test_body = std::make_unique<Body>(Body::load(state));
    // Save state
    auto saved_state = test_body->save();
    // Create new body from saved state
    Body loaded_body = Body::load(saved_state);
    // Verify loaded values
    const auto& X = loaded_body.getPosition();
    const auto& V = loaded_body.getVelocity();
    const auto& m = loaded_body.getMass();
    EXPECT_DOUBLE_EQ(X[0], 1.0);
    EXPECT_DOUBLE_EQ(X[1], 2.0);
    EXPECT_DOUBLE_EQ(X[2], 3.0);
    EXPECT_DOUBLE_EQ(V[0], 4.0);
    EXPECT_DOUBLE_EQ(V[1], 5.0);
    EXPECT_DOUBLE_EQ(V[2], 6.0);
    EXPECT_DOUBLE_EQ(m, 7.0);
 } 
--- a/tests/calc_test.cpp
+++ b/tests/calc_test.cpp
@ -0,0 +1,10 @@
 #include <gtest/gtest.h>
 #include "calc.hpp"
 #include <cmath>
 TEST(CalcTest, FormatSigFigs) {
    // Test positive numbers
    EXPECT_EQ(format_sig_figs(123.456, 3), "1.23e+02");
    EXPECT_EQ(format_sig_figs(123.456, 4), "1.235e+02");
    EXPECT_EQ(format_sig_figs(123.456, 5), "1.2346e+02");
 }
--- a/tests/simulator_test.cpp
+++ b/tests/simulator_test.cpp
@ -0,0 +1,54 @@
 // #include <gtest/gtest.h>
 // #include "simulator.hpp"
 // #include <vector>
 // #include <filesystem>
 // #include <memory>
 // class SimulatorTest : public testing::Test {
 // protected:
 //     SimulatorTest() : 
 //         bodies({
 //             Body(
 //                 Position{Decimal(0), Decimal(0), Decimal(0)},
 //                 Velocity{Decimal(0), Decimal(0), Decimal(0)},
 //                 Mass(1.0),
 //                 "body1"
 //             ),
 //             Body(
 //                 Position{Decimal(1), Decimal(0), Decimal(0)},
 //                 Velocity{Decimal(0), Decimal(1), Decimal(0)},
 //                 Mass(1.0),
 //                 "body2"
 //             )
 //         })
 //     {
 //         simulator = std::make_unique<Simulator>(
 //             bodies,
 //             Decimal(0.1),  // step_size
 //             1,            // steps_per_save
 //             std::filesystem::path("test_output.json"),  // output_file
 //             0,            // current_step
 //             true          // overwrite_output
 //         );
 //     }
 //     std::vector<Body> bodies;
 //     std::unique_ptr<Simulator> simulator;
 // };
 // TEST_F(SimulatorTest, Step) {
 //     // Take a step
 //     //TODO
 // }
 // TEST_F(SimulatorTest, TotalEnergy) {
 //     //TODO
 //     // Decimal initial_energy = simulator->total_energy();
 //     // simulator->step(Decimal(0.1));
 //     // Decimal new_energy = simulator->total_energy();
 //     // EXPECT_NEAR(initial_energy, new_energy, 1e-10);
 // }
--- a/tracjtory.out
+++ b/tracjtory.out
--- a/tracjtory.out2
+++ b/tracjtory.out2