Compare commits
11 Commits
7de52b0a3d
...
bd99aca317
| Author | SHA1 | Date | |
|---|---|---|---|
| bd99aca317 | |||
|
6cb8d8e450 | ||
|
|
b1b2e51507 | ||
|
|
be6229263d | ||
|
|
151e11076b | ||
|
|
63fe3a38b8 | ||
|
|
64bec5af20 | ||
|
|
e07bb7c4a1 | ||
|
|
16b71407d2 | ||
|
|
97cad94f6e | ||
|
|
b4d6062640 |
14
.gitignore
vendored
14
.gitignore
vendored
@ -4,4 +4,16 @@ __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
|
||||
@ -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()
|
||||
19
Cargo.toml
Normal file
19
Cargo.toml
Normal file
@ -0,0 +1,19 @@
|
||||
[package]
|
||||
name = "orbital_simulator"
|
||||
version = "0.1.0"
|
||||
edition = "2021"
|
||||
authors = ["thomas"]
|
||||
|
||||
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
|
||||
|
||||
[dependencies]
|
||||
bincode = "1.3"
|
||||
clap = { version = "4.5.39", features = ["derive"] }
|
||||
env_logger = "0.11.8"
|
||||
glam = { version = "0.30.4", features = ["serde"] }
|
||||
humantime = "2.2.0"
|
||||
indicatif = "0.17.11"
|
||||
log = "0.4.27"
|
||||
once_cell = "1.21.3"
|
||||
serde = { version = "1.0.219", features = ["derive"] }
|
||||
serde_json = "1.0.140"
|
||||
83
README.md
83
README.md
@ -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.
|
||||
165
animate.py
165
animate.py
@ -1,165 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.animation import FuncAnimation
|
||||
from collections import defaultdict
|
||||
import argparse
|
||||
import re
|
||||
|
||||
def parse_line(line):
|
||||
"""Parse a line from the simulation output."""
|
||||
# Skip comments and empty lines
|
||||
if line.startswith('#') or not line.strip():
|
||||
return None
|
||||
|
||||
# Parse the line format: body_name: X = x1m x2m x3m, V = v1m/s v2m/s v3m/s
|
||||
try:
|
||||
# Split on colon to separate name and data
|
||||
name_part, data_part = line.strip().split(':')
|
||||
name = name_part.strip()
|
||||
|
||||
# Split data part into position and velocity
|
||||
pos_part, vel_part = data_part.split(',')
|
||||
|
||||
# Extract position values
|
||||
pos_str = pos_part.split('=')[1].strip()
|
||||
pos_values = [float(x.replace('m', '').strip()) for x in pos_str.split() if x.strip()]
|
||||
|
||||
if len(pos_values) != 3:
|
||||
return None
|
||||
|
||||
return {
|
||||
'name': name,
|
||||
'position': tuple(pos_values)
|
||||
}
|
||||
except (ValueError, IndexError):
|
||||
return None
|
||||
|
||||
def read_output_file(filename):
|
||||
"""Read the simulation output file and organize data by body and time."""
|
||||
positions = defaultdict(list)
|
||||
times = [] # We'll use frame numbers as times since actual time isn't in output
|
||||
|
||||
with open(filename, 'r') as f:
|
||||
frame = 0
|
||||
for line in f:
|
||||
data = parse_line(line)
|
||||
if data is None:
|
||||
continue
|
||||
|
||||
name = data['name']
|
||||
pos = data['position']
|
||||
|
||||
# Store position
|
||||
positions[name].append((pos[0], pos[1]))
|
||||
frame += 1
|
||||
times.append(frame)
|
||||
|
||||
return positions, times
|
||||
|
||||
def center_positions(positions, center_body):
|
||||
"""Center all positions relative to the specified body."""
|
||||
if center_body not in positions:
|
||||
print(f"Warning: Center body '{center_body}' not found in simulation")
|
||||
return positions
|
||||
|
||||
centered_positions = defaultdict(list)
|
||||
for frame in range(len(next(iter(positions.values())))):
|
||||
# Get center body position for this frame
|
||||
center_pos = positions[center_body][frame]
|
||||
|
||||
# Shift all positions relative to center body
|
||||
for name, pos_list in positions.items():
|
||||
if frame < len(pos_list):
|
||||
pos = pos_list[frame]
|
||||
centered_pos = (pos[0] - center_pos[0], pos[1] - center_pos[1])
|
||||
centered_positions[name].append(centered_pos)
|
||||
|
||||
return centered_positions
|
||||
|
||||
def create_animation(positions, times, output_file=None, center_body=None):
|
||||
"""Create an animation of the bodies' orbits."""
|
||||
# Check if we have any data
|
||||
if not positions or not times:
|
||||
print("Error: No valid data found in the input file")
|
||||
return
|
||||
|
||||
# Center positions if requested
|
||||
if center_body:
|
||||
positions = center_positions(positions, center_body)
|
||||
|
||||
# Set up the figure and axis
|
||||
fig, ax = plt.subplots(figsize=(10, 10))
|
||||
|
||||
# Set equal aspect ratio
|
||||
ax.set_aspect('equal')
|
||||
|
||||
# Create scatter plots for each body
|
||||
scatters = {}
|
||||
for name in positions.keys():
|
||||
scatters[name], = ax.plot([], [], 'o-', label=name, alpha=0.7)
|
||||
|
||||
# Set up the plot
|
||||
ax.set_xlabel('X (m)')
|
||||
ax.set_ylabel('Y (m)')
|
||||
title = 'Orbital Simulation'
|
||||
if center_body:
|
||||
title += f' (centered on {center_body})'
|
||||
ax.set_title(title)
|
||||
ax.legend()
|
||||
|
||||
# Find the bounds of the plot
|
||||
all_x = []
|
||||
all_y = []
|
||||
for pos_list in positions.values():
|
||||
if pos_list: # Only process if we have positions
|
||||
x, y = zip(*pos_list)
|
||||
all_x.extend(x)
|
||||
all_y.extend(y)
|
||||
|
||||
if not all_x or not all_y:
|
||||
print("Error: No valid position data found")
|
||||
return
|
||||
|
||||
max_range = max(max(abs(min(all_x)), abs(max(all_x))),
|
||||
max(abs(min(all_y)), abs(max(all_y))))
|
||||
ax.set_xlim(-max_range, max_range)
|
||||
ax.set_ylim(-max_range, max_range)
|
||||
|
||||
def init():
|
||||
"""Initialize the animation."""
|
||||
for scatter in scatters.values():
|
||||
scatter.set_data([], [])
|
||||
return list(scatters.values())
|
||||
|
||||
def update(frame):
|
||||
"""Update the animation for each frame."""
|
||||
for name, scatter in scatters.items():
|
||||
if positions[name]: # Only update if we have positions
|
||||
x, y = zip(*positions[name][:frame+1])
|
||||
scatter.set_data(x, y)
|
||||
return list(scatters.values())
|
||||
|
||||
# Create the animation
|
||||
anim = FuncAnimation(fig, update, frames=len(times),
|
||||
init_func=init, blit=True,
|
||||
interval=50) # 50ms between frames
|
||||
|
||||
if output_file:
|
||||
anim.save(output_file, writer='ffmpeg', fps=30)
|
||||
else:
|
||||
plt.show()
|
||||
|
||||
def main():
|
||||
parser = argparse.ArgumentParser(description='Animate orbital simulation output')
|
||||
parser.add_argument('input_file', help='Input file from simulation')
|
||||
parser.add_argument('--output', '-o', help='Output video file (optional)')
|
||||
parser.add_argument('--center', '-c', help='Center the animation on this body')
|
||||
args = parser.parse_args()
|
||||
|
||||
positions, times = read_output_file(args.input_file)
|
||||
create_animation(positions, times, args.output, args.center)
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
22
config/jwst.json
Normal file
22
config/jwst.json
Normal file
@ -0,0 +1,22 @@
|
||||
{
|
||||
"bodies": [
|
||||
{
|
||||
"name": "Earth",
|
||||
"mass": 5.972e24,
|
||||
"position": [1.47095e11, 0, 0],
|
||||
"velocity": [0, 29290, 0]
|
||||
},
|
||||
{
|
||||
"name": "Sun",
|
||||
"mass": 1.989e30,
|
||||
"position": [0, 0, 0],
|
||||
"velocity": [0, 0, 0]
|
||||
},
|
||||
{
|
||||
"name": "JWST",
|
||||
"mass": 6500,
|
||||
"position": [149217067274.40607, 0, 0],
|
||||
"velocity": [0, 29729.784, 0]
|
||||
}
|
||||
]
|
||||
}
|
||||
78
src/body.cpp
78
src/body.cpp
@ -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
48
src/body.hpp
@ -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
147
src/calc.cpp
@ -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
22
src/calc.hpp
@ -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
|
||||
19
src/config.rs
Normal file
19
src/config.rs
Normal file
@ -0,0 +1,19 @@
|
||||
use crate::types;
|
||||
use serde::{Deserialize, Serialize};
|
||||
|
||||
|
||||
#[derive(Debug, Deserialize, Serialize)]
|
||||
pub struct Body {
|
||||
pub name: String,
|
||||
pub mass: types::Mass,
|
||||
pub position: types::Position,
|
||||
pub velocity: types::Velocity,
|
||||
|
||||
#[serde(default)]
|
||||
pub acceleration: types::Acceleration,
|
||||
}
|
||||
|
||||
#[derive(Debug, Deserialize)]
|
||||
pub struct ConfigFile {
|
||||
pub bodies: Vec<Body>,
|
||||
}
|
||||
167
src/main.cpp
167
src/main.cpp
@ -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;
|
||||
}
|
||||
}
|
||||
122
src/main.rs
Normal file
122
src/main.rs
Normal file
@ -0,0 +1,122 @@
|
||||
// Standard library
|
||||
use std::error::Error;
|
||||
use std::fs::File;
|
||||
use std::io::BufReader;
|
||||
use std::path::Path;
|
||||
use std::time::{Duration,Instant};
|
||||
|
||||
// External crates
|
||||
use clap::Parser;
|
||||
use log::{info, debug};
|
||||
use indicatif::{ProgressBar, ProgressStyle};
|
||||
use humantime;
|
||||
|
||||
// Internal modules
|
||||
mod types;
|
||||
mod config;
|
||||
mod simulator;
|
||||
|
||||
// Specific uses from modules
|
||||
use crate::simulator::Simulation;
|
||||
use crate::types::{norm_time, real_body};
|
||||
|
||||
|
||||
#[derive(Parser, Debug)]
|
||||
#[command(
|
||||
version,
|
||||
about="Orbital mechanics simulator",
|
||||
long_about = "Given initial conditions you provide to --config, \
|
||||
this program will numerically integrate and determinate their \
|
||||
paths based off Newton's law of gravity.")]
|
||||
struct Args {
|
||||
///Config file for initial conditions
|
||||
#[arg(short, long)]
|
||||
config: String,
|
||||
|
||||
/// Time to run simulation for (e.g. 10s, 5m, 2h, 100d)
|
||||
#[arg(short, long, value_parser = humantime::parse_duration)]
|
||||
time: Duration,
|
||||
|
||||
///Step size for simulation (seconds)
|
||||
#[arg(short, long, default_value_t = 10.0)]
|
||||
step_size: f64,
|
||||
|
||||
///How often to update progress bar and save (seconds)
|
||||
#[arg(short = 'P', long, default_value_t = 1000)]
|
||||
steps_per_save: usize,
|
||||
|
||||
/// Filename to save trajectory to
|
||||
#[arg(short, long)]
|
||||
output_file: String,
|
||||
}
|
||||
|
||||
fn read_config<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);
|
||||
})
|
||||
);
|
||||
}
|
||||
@ -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();
|
||||
}
|
||||
@ -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;
|
||||
};
|
||||
113
src/simulator.rs
Normal file
113
src/simulator.rs
Normal file
@ -0,0 +1,113 @@
|
||||
use std::fs::OpenOptions;
|
||||
use std::io::BufWriter;
|
||||
|
||||
use serde::Serialize;
|
||||
use log::{debug, trace};
|
||||
use glam::DVec3;
|
||||
|
||||
use crate::config::Body;
|
||||
use crate::types::{norm_mass, norm_pos, norm_vel, norm_time, real_pos, real_vel, real_time, real_body};
|
||||
|
||||
pub struct Simulation {
|
||||
pub bodies: Vec<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));
|
||||
}
|
||||
}
|
||||
}
|
||||
105
src/types.rs
Normal file
105
src/types.rs
Normal file
@ -0,0 +1,105 @@
|
||||
#![allow(unused)]
|
||||
use glam::DVec3;
|
||||
use once_cell::sync::Lazy;
|
||||
|
||||
use crate::config::Body;
|
||||
|
||||
pub type Position = DVec3;
|
||||
pub type Velocity = DVec3;
|
||||
pub type Acceleration = DVec3;
|
||||
pub type Mass = f64;
|
||||
pub type Time = f64;
|
||||
|
||||
|
||||
// Constants
|
||||
const EARTH_RADIUS: f64 = 6.378e6; //meters
|
||||
const EARTH_MASS: f64 = 5.972e24; //kg
|
||||
const EARTH_ORBITAL_VELOCITY: f64 = 2.9780e4; //m/s
|
||||
const AU: f64 = 1.49597870700e11; //meters
|
||||
|
||||
const MOON_MASS: Mass = 7.34767309e22; // kg
|
||||
const MOON_ORBITAL_VELOCITY: f64 = 1.022e3; //m/s relative to earth
|
||||
|
||||
const SUN_MASS: f64 = 1.989e30; //kg
|
||||
const SUN_RADIUS: f64 = 6.957e8; //meters
|
||||
|
||||
const G: f64 = 6.67430e-11;
|
||||
const R_0: f64 = EARTH_RADIUS;
|
||||
const M_0: f64 = EARTH_MASS;
|
||||
|
||||
static T_0: Lazy<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,
|
||||
}
|
||||
}
|
||||
@ -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
44
test.py
@ -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()
|
||||
@ -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);
|
||||
}
|
||||
@ -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");
|
||||
}
|
||||
@ -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);
|
||||
// }
|
||||
|
||||
BIN
tracjtory.out
Normal file
BIN
tracjtory.out
Normal file
Binary file not shown.
BIN
tracjtory.out2
Normal file
BIN
tracjtory.out2
Normal file
Binary file not shown.
Loading…
x
Reference in New Issue
Block a user