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

Trying a C++ rewrite

Browse Source
This commit is contained in:
Thomas Faour 2025-06-01 19:38:37 -04:00
parent 58eea0f264
commit 4b21fb75e8
11 changed files with 669 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitignore vendored
View File

@ -1,3 +1,7 @@
last_checkpoint.npz last_checkpoint.npz
*.pyc *.pyc
__pycache__ __pycache__
/build
/cmake-build-*
CMakeFiles/
*.cmake

47
CMakeLists.txt Normal file
View File

@ -0,0 +1,47 @@
cmake_minimum_required(VERSION 3.10)
project(orbital_simulator)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
# Find required packages
find_package(Boost REQUIRED COMPONENTS system)
# Add source files
set(SOURCES
src/main.cpp
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 executable
add_executable(orbital_simulator ${SOURCES} ${HEADERS})
# Link libraries
target_link_libraries(orbital_simulator PRIVATE
Boost::system
${Boost_LIBRARIES}
)
# Include directories
target_include_directories(orbital_simulator PRIVATE
${CMAKE_SOURCE_DIR}
${Boost_INCLUDE_DIRS}
)
# 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()

84
README.md
View File

@ -1 +1,83 @@
Of course the difficult thing about an orbital simulator is the magnitude problem. Such big differences in magnitude! I'll probably need to calculate forces, etc, within a specific body. And that body can then use units that are proportional to only it. # 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.

66
src/body.cpp Normal file
View File

@ -0,0 +1,66 @@
#include "body.hpp"
#include "calc.hpp"
#include <cmath>
#include <boost/multiprecision/cpp_dec_float.hpp>
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)};
}
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 boost::multiprecision::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 boost::multiprecision::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]), 3) + "m ";
vel_str += format_sig_figs(real_vel(V[i]), 3) + "m/s ";
}
return name + ": X = " + pos_str + ", V = " + vel_str;
}

44
src/body.hpp Normal file
View File

@ -0,0 +1,44 @@
#pragma once
#include "units.hpp"
#include <string>
#include <tuple>
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; }
// String representation
std::string toString() const;
private:
Position X;
Velocity V;
Acceleration A;
Mass m;
std::string name;
Decimal _speed() const; // Internal speed calculation
};

131
src/calc.cpp Normal file
View File

@ -0,0 +1,131 @@
#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 = boost::multiprecision::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) {
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();
int steps_per_second = elapsed > 0 ? iteration / elapsed : 0;
std::cout << "\r" << bar << " " << std::fixed << std::setprecision(2)
<< percent << "% " << steps_per_second << " steps/s" << 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 Normal file
View File

@ -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 = 50);
// 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

43
src/main.cpp Normal file
View File

@ -0,0 +1,43 @@
#include "simulator.hpp"
#include "units.hpp"
#include <iostream>
int main() {
try {
// Create bodies
std::vector<Body> bodies;
// Earth at origin
Position earth_pos{Decimal(0), Decimal(0), Decimal(0)};
Velocity earth_vel{Decimal(0), Decimal(0), Decimal(0)};
bodies.emplace_back(earth_pos, earth_vel, EARTH_MASS, "Earth");
// Moon at average distance (384,400 km) with orbital velocity
// Using average orbital velocity of 1.022 km/s
Position moon_pos{Decimal("384400000"), Decimal(0), Decimal(0)}; // 384,400 km in meters
Velocity moon_vel{Decimal(0), Decimal("1022"), Decimal(0)}; // 1.022 km/s in m/s
bodies.emplace_back(moon_pos, moon_vel, MOON_MASS, "Moon");
// Create simulator
// Step size: 100 seconds (small enough for accuracy)
// Steps per save: 1000 (save every ~27.8 hours)
// Total steps: 27.3 days * 24 hours * 3600 seconds / 100 seconds per step
Decimal step_size = Decimal("1"); // 100 seconds per step
int steps_per_save = 1000; // Save every 1000 steps
int total_steps = 27.3 * 24 * 3600 / 1; // Steps for one orbit
Simulator sim(bodies, step_size, steps_per_save, std::filesystem::path("output.txt"));
// Run simulation
std::cout << "Starting simulation of one lunar orbit...\n";
std::cout << "Step size: " << step_size << " seconds\n";
std::cout << "Total steps: " << total_steps << " (approximately 27.3 days)\n";
sim.run(total_steps);
std::cout << "\nSimulation complete!\n";
return 0;
} catch (const std::exception& e) {
std::cerr << "Error: " << e.what() << std::endl;
return 1;
}
}

149
src/simulator.cpp Normal file
View File

@ -0,0 +1,149 @@
#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();
}
// Write initial header with masses
out = std::ofstream(output_file, std::ios::app);
if (!out) {
throw std::runtime_error("Failed to open output file: " + output_file.string());
}
for (const auto& body : bodies) {
out << body.getName() << ": " << body.getMass() << "\n";
}
out << "\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();
current_step++;
if (current_step % steps_per_save == 0) {
print_progress_bar(i, steps, start_time);
checkpoint();
}
}
}
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] = bodies[i].getPosition()[k] - bodies[j].getPosition()[k];
}
// Calculate distance
Decimal dist = boost::multiprecision::sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
// Calculate force magnitude using Newton's law of gravitation
// F = G * m1 * m2 / r^2
Decimal force_magnitude = G * bodies[i].getMass() * bodies[j].getMass() / (dist * dist);
// Calculate acceleration for both bodies
// a = F/m
Decimal acc_magnitude_i = force_magnitude / bodies[i].getMass();
Decimal acc_magnitude_j = force_magnitude / bodies[j].getMass();
// Convert to vector form
Acceleration acc_i, acc_j;
for (int k = 0; k < 3; ++k) {
acc_i[k] = -vec[k] * acc_magnitude_i / dist; // Negative because force is attractive
acc_j[k] = vec[k] * acc_magnitude_j / dist; // Positive because force is attractive
}
// Add to current accelerations
Acceleration current_acc_i = bodies[i].getAcceleration();
Acceleration current_acc_j = bodies[j].getAcceleration();
for (int k = 0; k < 3; ++k) {
current_acc_i[k] += acc_i[k];
current_acc_j[k] += acc_j[k];
}
bodies[i].setAcceleration(current_acc_i);
bodies[j].setAcceleration(current_acc_j);
}
}
}
void Simulator::move_bodies() {
for (auto& body : bodies) {
body.step(step_size);
}
}
void Simulator::checkpoint() {
// Open file in append mode
std::ofstream out(output_file, std::ios::app);
if (!out) {
throw std::runtime_error("Failed to open output file: " + output_file.string());
}
// Write timestamp and step number
out << "Time: " << (current_step * step_size) << "s\n";
// Write positions and velocities for each body
for (const auto& body : bodies) {
out << body.getName() << " : ";
out << " ";
for (int i = 0; i < 3; ++i) {
out << real_pos(body.getPosition()[i]) << "m ";
}
out << " : ";
for (int i = 0; i < 3; ++i) {
out << real_vel(body.getVelocity()[i]) << "m/s ";
}
out << "\n";
}
}

35
src/simulator.hpp Normal file
View File

@ -0,0 +1,35 @@
#pragma once
#include "body.hpp"
#include <vector>
#include <string>
#include <filesystem>
#include <memory>
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);
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;
};

44
src/units.hpp Normal file
View File

@ -0,0 +1,44 @@
#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 = Decimal("5972e21"); // kg
const Decimal EARTH_RADIUS = Decimal("6378e3"); // meters
const Decimal EARTH_ORBITAL_VELOCITY = Decimal("29780"); // m/s
const Decimal AU = Decimal("149597870700"); // meters
const Decimal MOON_MASS = Decimal("734767309e14");
const Decimal MOON_ORBITAL_VELOCITY = Decimal("1022"); // m/s relative to earth
const Decimal SUN_MASS = Decimal("1989e27"); // kg
const Decimal SUN_RADIUS = Decimal("6957e5"); // meters
const Decimal PI = Decimal("3.14159265358979323846264338327950288419716939937510");
// Normalizing constants
const Decimal G = Decimal("6.67430e-11");
const Decimal r_0 = EARTH_RADIUS;
const Decimal m_0 = Decimal("5.972e24");
const Decimal t_0 = boost::multiprecision::sqrt((r_0 * r_0 * r_0) / (G * m_0));
// 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); }