tfaour/OrbitalSimulator
1
0
Fork 0
Code Issues 4 Pull Requests Actions Packages Projects Releases Wiki Activity

initial rust

Browse Source
This commit is contained in:
Thomas Faour 2025-06-06 20:04:43 -04:00
parent b4d6062640
commit 97cad94f6e
19 changed files with 41 additions and 1418 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitignore vendored
View File

@ -11,3 +11,9 @@ rust/Cargo.lock
out.out out.out
output.json output.json
*.png *.png
# Added by cargo
/target
Cargo.lock

85
CMakeLists.txt
View File

@ -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()

10
Cargo.toml Normal file
View File

@ -0,0 +1,10 @@
[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]
clap = { version = "4.5.39", features = ["derive"] }

83
README.md
View File

@ -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<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.

201
animate.py
View File

@ -1,201 +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 calculate_l2_point(sun_pos, earth_pos):
"""Calculate L2 Lagrange point position."""
# Calculate distance between bodies
dx = earth_pos[0] - sun_pos[0]
dy = earth_pos[1] - sun_pos[1]
r = np.sqrt(dx*dx + dy*dy)
# Calculate mass ratio (using approximate values)
mu = 5.972e24 / (1.989e30 + 5.972e24) # Earth mass / (Sun mass + Earth mass)
# Calculate L2 position (beyond Earth)
L2_x = sun_pos[0] + (1 + mu**(1/3)) * dx
L2_y = sun_pos[1] + (1 + mu**(1/3)) * dy
return (L2_x, L2_y)
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 (L2 centered)'
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
# Calculate initial L2 point for reference
if 'Sun' in positions and 'Earth' in positions:
initial_l2 = calculate_l2_point(positions['Sun'][0], positions['Earth'][0])
# Set a fixed zoom level around L2 (adjust the factor to change zoom)
zoom_factor = 2e9 # 2 million km view
ax.set_xlim(-zoom_factor, zoom_factor)
ax.set_ylim(-zoom_factor, zoom_factor)
# Create L2 point scatter plot at origin
l2_scatter, = ax.plot([0], [0], 'gx', markersize=10, label='L2')
scatters['L2'] = l2_scatter
else:
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."""
if 'Sun' in positions and 'Earth' in positions:
# Calculate L2 point for current frame
l2_point = calculate_l2_point(positions['Sun'][frame], positions['Earth'][frame])
# Update positions relative to L2
for name, scatter in scatters.items():
if name == 'L2':
scatter.set_data([0], [0]) # Keep L2 at origin
elif positions[name]: # Only update if we have positions
x, y = zip(*positions[name][:frame+1])
# Shift positions relative to L2
x = [pos - l2_point[0] for pos in x]
y = [pos - l2_point[1] for pos in y]
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()

158
python/plot_potential.py
View File

@ -1,158 +0,0 @@
import json
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.colors import LogNorm
# Constants from our Rust implementation
G = 6.67430e-11 # Gravitational constant
R_0 = 149597870700 # Earth radius (normalization length)
M_0 = 5.972e24 # Earth mass (normalization mass)
T_0 = 5060.0 # Characteristic time
def load_bodies(config_file):
"""Load bodies from config file."""
with open(config_file, 'r') as f:
data = json.load(f)
return data['bodies']
def calculate_potential(x, y, bodies):
"""Calculate gravitational potential at point (x,y)."""
potential = 0.0
for body in bodies:
dx = x - body['position'][0]
dy = y - body['position'][1]
r = np.sqrt(dx*dx + dy*dy)
if r > 0: # Avoid division by zero
potential -= G * body['mass'] / r
return potential
def calculate_lagrange_points(bodies):
"""Calculate approximate positions of Lagrange points."""
if len(bodies) != 2:
return []
# Find primary and secondary bodies
if bodies[0]['mass'] > bodies[1]['mass']:
primary, secondary = bodies[0], bodies[1]
else:
primary, secondary = bodies[1], bodies[0]
# Calculate distance between bodies
dx = secondary['position'][0] - primary['position'][0]
dy = secondary['position'][1] - primary['position'][1]
r = np.sqrt(dx*dx + dy*dy)
# Calculate mass ratio
mu = secondary['mass'] / (primary['mass'] + secondary['mass'])
# Calculate Lagrange points
L1 = np.array([primary['position'][0] + (1 - mu**(1/3)) * dx,
primary['position'][1] + (1 - mu**(1/3)) * dy])
L2 = np.array([primary['position'][0] + (1 + mu**(1/3)) * dx,
primary['position'][1] + (1 + mu**(1/3)) * dy])
L3 = np.array([primary['position'][0] - (1 + 5*mu/12) * dx,
primary['position'][1] - (1 + 5*mu/12) * dy])
L4 = np.array([primary['position'][0] + 0.5 * dx - 0.866 * dy,
primary['position'][1] + 0.5 * dy + 0.866 * dx])
L5 = np.array([primary['position'][0] + 0.5 * dx + 0.866 * dy,
primary['position'][1] + 0.5 * dy - 0.866 * dx])
return [L1, L2, L3, L4, L5]
def plot_potential(config_file, grid_size=200, extent=1.5):
"""Plot gravitational potential field."""
# Load bodies
bodies = load_bodies(config_file)
# Find Earth and Sun
earth = next((body for body in bodies if body['name'] == 'Earth'), None)
sun = next((body for body in bodies if body['name'] == 'Sun'), None)
if not (earth and sun):
print("Error: Earth and Sun must be present in the config file")
return
# Calculate Lagrange points for Earth-Sun system
lagrange_points = calculate_lagrange_points([sun, earth])
# Create coordinate grid centered on Earth
x = np.linspace(earth['position'][0] - extent*R_0, earth['position'][0] + extent*R_0, grid_size)
y = np.linspace(earth['position'][1] - extent*R_0, earth['position'][1] + extent*R_0, grid_size)
X, Y = np.meshgrid(x, y)
# Calculate potential at each grid point
Z = np.zeros_like(X)
for i in range(grid_size):
for j in range(grid_size):
Z[i,j] = calculate_potential(X[i,j], Y[i,j], bodies)
# Take absolute value for logarithmic scale
Z_abs = np.abs(Z)
# Calculate reference potential at Earth's position
earth_potential = abs(calculate_potential(earth['position'][0], earth['position'][1], [sun]))
# Set scale to focus on Earth-Sun system
vmin = earth_potential / 10 # Show potential variations within an order of magnitude
vmax = earth_potential * 10
# Create plot
plt.figure(figsize=(12, 10))
ax = plt.gca()
# Plot potential field with logarithmic scale
im = ax.imshow(Z_abs, extent=[x[0], x[-1], y[0], y[-1]],
origin='lower', cmap='viridis', norm=LogNorm(vmin=vmin, vmax=vmax))
# Add colorbar
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
plt.colorbar(im, cax=cax, label='|Gravitational Potential| (J/kg)')
# Plot bodies with different colors
body_colors = {
'Sun': 'yellow',
'Earth': 'blue',
'Moon': 'gray',
'Mars': 'red',
'Venus': 'orange',
'Mercury': 'brown',
'Jupiter': 'orange',
'Saturn': 'gold',
'Uranus': 'lightblue',
'Neptune': 'blue'
}
for body in bodies:
color = body_colors.get(body['name'], 'white') # Default to white if name not found
ax.plot(body['position'][0], body['position'][1], 'o',
color=color, markersize=15, label=body['name'])
# Plot Lagrange points
for i, point in enumerate(lagrange_points, 1):
ax.plot(point[0], point[1], 'r+', markersize=12, label=f'L{i}')
# Customize plot
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_title('Gravitational Potential Field with Lagrange Points (Log Scale)')
ax.legend()
# Save plot
plt.savefig('potential_field.png', dpi=300, bbox_inches='tight')
plt.show()
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description='Plot gravitational potential field from config file')
parser.add_argument('config_file', help='Path to config file')
parser.add_argument('--grid-size', type=int, default=200, help='Grid size for potential calculation')
parser.add_argument('--extent', type=float, default=1.5, help='Plot extent in Earth radii')
args = parser.parse_args()
plot_potential(args.config_file, args.grid_size, args.extent)

78
src/body.cpp
View File

@ -1,78 +0,0 @@
#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;
}

48
src/body.hpp
View File

@ -1,48 +0,0 @@
#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
};

147
src/calc.cpp
View File

@ -1,147 +0,0 @@
#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

22
src/calc.hpp
View File

@ -1,22 +0,0 @@
#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

167
src/main.cpp
View File

@ -1,167 +0,0 @@
#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;
}
}

25
src/main.rs Normal file
View File

@ -0,0 +1,25 @@
use clap::Parser;
#[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,
///Step size for simulation (seconds)
#[arg(short, long, default_value_t = 10)]
step_size: u8,
}
fn main() {
let args = Args::parse();
println!("Loading initial parameters from {}", args.config);
}

121
src/simulator.cpp
View File

@ -1,121 +0,0 @@
#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();
}

43
src/simulator.hpp
View File

@ -1,43 +0,0 @@
#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;
};

46
src/units.hpp
View File

@ -1,46 +0,0 @@
#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)); }

44
test.py
View File

@ -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()

111
tests/body_test.cpp
View File

@ -1,111 +0,0 @@
#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);
}

10
tests/calc_test.cpp
View File

@ -1,10 +0,0 @@
#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");
}

54
tests/simulator_test.cpp
View File

@ -1,54 +0,0 @@
// #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);
// }