From f52ecf4889468aeb97d00b2a06efe8392faec4c6 Mon Sep 17 00:00:00 2001 From: Thomas Faour Date: Mon, 2 Jun 2025 08:27:45 -0400 Subject: [PATCH] Added config file for planets --- CMakeLists.txt | 16 +- animate.py | 135 +++++++++++++ config/planets.json | 22 +++ config/planets_fake.json | 16 ++ .../orbits/__pycache__/body.cpython-312.pyc | Bin 1644 -> 3915 bytes .../__pycache__/simulator.cpython-312.pyc | Bin 4035 -> 6147 bytes src/body.cpp | 4 +- src/calc.cpp | 26 ++- src/calc.hpp | 2 +- src/main.cpp | 179 +++++++++++++++--- src/simulator.cpp | 51 ++--- src/simulator.hpp | 1 + src/units.hpp | 24 +-- 13 files changed, 386 insertions(+), 90 deletions(-) create mode 100644 animate.py create mode 100644 config/planets.json create mode 100644 config/planets_fake.json diff --git a/CMakeLists.txt b/CMakeLists.txt index e18eabe..8dc063e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,8 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) # Find required packages -find_package(Boost REQUIRED COMPONENTS system) +find_package(Boost REQUIRED COMPONENTS program_options) +find_package(nlohmann_json REQUIRED) # Add source files set(SOURCES @@ -27,15 +28,16 @@ set(HEADERS add_executable(orbital_simulator ${SOURCES} ${HEADERS}) # Link libraries -target_link_libraries(orbital_simulator PRIVATE - Boost::system - ${Boost_LIBRARIES} +target_link_libraries(orbital_simulator + PRIVATE + Boost::program_options + nlohmann_json::nlohmann_json ) # Include directories -target_include_directories(orbital_simulator PRIVATE - ${CMAKE_SOURCE_DIR} - ${Boost_INCLUDE_DIRS} +target_include_directories(orbital_simulator + PRIVATE + ${CMAKE_SOURCE_DIR}/src ) # Optional: Enable ncurses for terminal plotting diff --git a/animate.py b/animate.py new file mode 100644 index 0000000..2ae5a3d --- /dev/null +++ b/animate.py @@ -0,0 +1,135 @@ +#!/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'] + positions[name].append((pos[0], pos[1])) + frame += 1 + times.append(frame) + + return positions, times + +def create_animation(positions, times, output_file=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 + + # 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)') + ax.set_title('Orbital Simulation') + 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)') + args = parser.parse_args() + + positions, times = read_output_file(args.input_file) + create_animation(positions, times, args.output) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/config/planets.json b/config/planets.json new file mode 100644 index 0000000..6e44305 --- /dev/null +++ b/config/planets.json @@ -0,0 +1,22 @@ +{ + "planets": [ + { + "name": "Earth", + "mass": 5.972e24, + "position": [149597870700, 0, 0], + "velocity": [0, 29780, 0] + }, + { + "name": "Moon", + "mass": 7.34767309e22, + "position": [149982270700, 0, 0], + "velocity": [0, 30802, 0] + }, + { + "name": "Sun", + "mass": 1.989e30, + "position": [0, 0, 0], + "velocity": [0, 0, 0] + } + ] +} \ No newline at end of file diff --git a/config/planets_fake.json b/config/planets_fake.json new file mode 100644 index 0000000..d16adeb --- /dev/null +++ b/config/planets_fake.json @@ -0,0 +1,16 @@ +{ + "planets": [ + { + "name": "Earth", + "mass": 7.342e24, + "position": [0, 0, 0], + "velocity": [0, -1022, 0] + }, + { + "name": "Moon", + "mass": 7.34767309e24, + "position": [384400000, 0, 0], + "velocity": [0, 1022, 0] + } + ] +} diff --git a/orbiter/orbits/__pycache__/body.cpython-312.pyc b/orbiter/orbits/__pycache__/body.cpython-312.pyc index 0ab85630d5607ded55a3ea3011f3a7d654c226b3..9991296f41bc9a4de3b319d5c022906addf00a77 100644 GIT binary patch literal 3915 zcmd59T}&Lud1n7ue$O0WgY7hHhXCs&#;#+G6W0z1yGdw}BUf&8RkgX?EO7Xah1~_r z`H=j;18x$5)kxqtm2n!$L8K~=ezc$R(3iufIM&vYQnfFWd&_}w^V*U2 zn{Vd(X6BpkclPhq)fEKVj~|{Iydo3wFB}w`6d;9eK;{V}jHXDM#wbNxNC`2KiA+jK ze-&UxW=cvHozNWW#L zMF2+GE`W8c7GV8kRIm4> zZi`K3^a>4Hcib_xL@Mq$8VIh%hbJ}5{wV3jQ`$$s#d$pXfk!WTH0?>l@wAx$9fu9| zE5cX!5Wqai(kw9m^IKMc5$c6EH(WOWR``-G2NdAJ9eB%%Fe=)NWF;7dt3dl=<_MPv zCBLk-jaW|7O;)S&5`OscRwIl4zriHbI}R*yt~IB>=CgE0Hi$rX9z9 zJlbVfA*mXGE_;$=rusLTH)#$Eg+gdQoo*kr(q_BcAGgNr_7U3}G83-T9w-sqAU&3f zyO!l z7k9k?@>70?i6SOW5MAlZ2x)B%6`GOK?T(gF-q1eKeg*xj+C}Jf)su#-Wa&V19k(zX2NUHaX=Hy-Tm&DZ>9Mtlsx&E}KWlXJgcR2H2@^Xrb4 zmsT2bbtm$X&Yar$bVDXh@6xB<-d+do^7;Pm-Rk|_`@8q4JLlLQVB$G;2M=5Z8aD96 zAWQjU9~P+@7+xSYcq67O;gabIBc{z067Cr!EP6Jha1EVpfk_8bALu4hH#>1{;^WDM z*KeJ;apLC5)kt$rZGJXj@OU0DrvZ6hNvQleK|vv*pu8}8l`QGFM!#*P4Qso|_7ONM zQ^tiJF?AR?m=3~1Ff%dde!r)6bx&)4kG>kwbE^KV;OJ~za2Jqy;)YuO+31FqQxvO4 z2BHz7w|r~#4ux+F*f?lKijsz*vEi9vy0!m zL_Jv73x072CtPg4({j7zVV%D4{#@7G=%0EX)ag)TzK^u8M%wd{qdE2HQ}* zgYfY9>41)>JNpDq1U2v#z)NeOjpZ&=>1Q*p9`pyh6TwLUThd7A-oFJ(sU_VoVQdx( zupB;ZwzXn%2iphr{3a>RfenSeH8eNWb$mi0L9#HI8WvRXDYb#9LGyL8UVCj z^fT+L!LQxFU0J%4Z+Y`>SAO^DT*H}sq${U(al`e&e-Y6e@P&)N4Tj-W7)IJ+V=2U| z3}bXGo(d?C`I(h~Y4v>YruWo~!DHT&-h5+e28*6dF`?B*Ng+yiZcpenmiE%S$(3Ya??WFI)QHS zGVfM;cy~`M;N(UR^({z0PuCJjH31OJhZ&~Cax5c0$S3$RG$h5!Hn delta 928 zcmZ8f&rcIU6n-;1yWRdmOG6b0X*SZF8&{w7~{pPdR7t-g_@9f^1UsH_>%qJH#6V7`QGg84g5@X?WNNRz*xU^ zcVUfhc3osMo7b;Q`lf68mTUPjH|E=}&49x!NV5sj+U2wXz>U|r6Wb^1uuwdxEc9XPKr0s;cY!sYUk~T+?6qu-LHL zp!4xwJ#5V1J}noOVd^{y8A4uV5V36GxwjH0_cw`V^;=_zW%UQ+5#MKeKn&zMx^lH3 zk99z(wC zhC)lA=mDr0K{r7TAxsv_g&^>~+Co*TOEhDau+5sfLXutrsubOsqRWneDh2e6nsC6m zX}7GHncW^~0j4)EqV#U2O{2qZW{xyk7Fd}B&djm4fviQWWpQ(G+xrEWI%iOMClLw{ zAGM-;4NYWA^^jMVLUXyYSPTAdgx*TW)>*5+?R$PzE|Tz UK=*f;_yiM2BF6>~0UbO^`N;N>rLA85*Y)c`nVx zd3f__AtT1cj4f`n^g`O6vB&Mei>fW{$T;Ipi?^qxcR1i2sy*#Ojw3G;)p?t!(hi6F z5fYczxR~n+vN383R}=Z0b1}~Rj&dctnn~Y3<_Ws4l!Ku~cyAJo#HmW+oJ!-olQ@amRe5LZS@{MB2N8B)368j1 zm2Q)`M|A=1Qe~iu>IUjv^Td=t9d7W?{129Gl00r_jiRt~#23lJ3GJe2c*S13w`kM5 z@^)iK=$odnLO8bLKoXB79wghbB!Xl=mITXJbimuzsC-fQE%`my-kFjjS@aq6V#Xyg z9#%00i1C3KR@^9gkaQvO8vhgr!%r!3F+=C9+c;}ueZ#=&KLv7&eBMnWgH_jV?>E+k zhmKI)6*T@}8!{f+F7b~Y@_PZhN^N*hCK+W%;Hm@In9JmoY27-4DX58jg7xEy0VIP+ ztm9bw;h9VxGR*Kux(mo)XRuSo-yP@Z^m~41iJ$L-$(Zo8uPA8eViFsIw<#tVORQPG zO!@qG3E}_YsO{wf3GHw5D9uhFe4bx!!LlG z%z7IzB`%Ni)+o|;K4!^}5`>f>r;xdMF(LqCM`pZ8enJ;9N>2@y79avi?gT!_qQjVX z1Og(TOP8ND5?)#P~U882}_Th$h>=hT#@BPhx(j?(g z2!St#IWQ%x!UY|&oi}MpN1S6O83x$o*-EzHzLv`1Ryj5fVg~d&mgX(8A2?;awl3{^ zgLhYqNRx=pi|aFWS=p1LHTi1vp*(B}k(wNNDEC`@-$!y^)mxK?9?5-m#lNRSYf7{_ zyEigbgGZTqe_~g;w66?RC#u^0_#@@QfeXB{Uf9q0#2J)UAr_lhHN@zXP8**}z6%v@ zLlH=cE75$XERRd9#OIF-GVcU8*03CSMc>3x;hzA&X8-`9V^b8qYIUYQ1?#TDI51vd@3gGG>aA?;@^!s!d?=8Oe;`7W@J&cU z;NKuIZte-f7hsC0W@A9iU^i21 z7CQSj@^le;jg~a&4IN7m;fce7D4PULA>4szXU|q)yaDri7i8=fc|0)u)|=aJe!+7A zsovAOK6~f#uP!%0Dog9m1_ABoHsU=~UE5VA)7M(<_!Mkj0Vj7-*)$w$iI*IxTqW$?far##@{68~;G;j%j7ZkG{sumP)s(YYAkYU_fY z*OslG9^*Q>*0F-xnVbT?5a?v`L__-pOzEeA+#+>Br7+j!Y?zKqZ@#T)5dB5(505=~yuGgG*p%cWTu=UmBj zLf11|{$@^XeT|`(;Uj}&8VN#;;VB;wYljv9u?86u{OLafa+@4*0nyR8Ku(X=W8(*1 z0kQXSITYD?Wm{^9AT(?w5^dO#bCBpr!-wK rA&&G>APt^U`Y{RpjdXuP&VNh-;m=5LpZNAk;9rg+rK3%P)N=HHA>j6x;E*Q!I zQiJ%?4Xaiv9FUH6wQARwFm!zM4rwiLvJFKri4CS;6V((o(7=v){bW3=E(^BSly*AV zX3UN>!PM(fL*IQ!UWk!jVR=J(K)f8y-4{EFTzx}&!*#zEcTfXA&$y=ncROkU8iZTj z{D#0Ssh!DVBDZMQqTG_ZR=ZbL>21g?UhDDp=n(1S-80^AdTUbi<>I5VTXX%M`La`- z3;7&lRmP(H3Z^1x2>9+S=55GlGbfy3VhqE&_eo~GL!7!#I6Q?Q0kKowyfUI@$F`)e zYI@%*DYGjjFhy%gv(*|a60wd~+OO6{h_a)n^S?n}-vkXEOArupKLq)#0 zZaGZ#EGyT-tr7ZHPD1JZMNKiwC*4U_{s%ETtA>Ry`ciCksBg#1w4 z7atok>-TE;Iw!Ltwb0c_Z#kr!H>d3_e|Us#eVx%PTd8EAyqr zYSm!_NTW-R<+_z(tyHyHo_kn=!$)8#%`o diff --git a/src/body.cpp b/src/body.cpp index 5781641..9d5dc9a 100644 --- a/src/body.cpp +++ b/src/body.cpp @@ -37,7 +37,7 @@ Decimal Body::dist_from_o() const { for (const auto& x : X) { sum += x * x; } - return boost::multiprecision::sqrt(sum); + return std::sqrt(sum); } Decimal Body::ke() const { @@ -49,7 +49,7 @@ Decimal Body::_speed() const { for (const auto& v : V) { sum += v * v; } - return boost::multiprecision::sqrt(sum); + return std::sqrt(sum); } std::string Body::speed() const { diff --git a/src/calc.cpp b/src/calc.cpp index 37b305b..7526c65 100644 --- a/src/calc.cpp +++ b/src/calc.cpp @@ -16,7 +16,7 @@ std::vector> calculate_distances(const std::vector start_time, int length) { +void print_progress_bar(int iteration, int total, std::chrono::time_point start_time, int length, Decimal step_size) { float percent = (float)iteration / total * 100; int filled_length = length * iteration / total; @@ -47,10 +47,26 @@ void print_progress_bar(int iteration, int total, std::chrono::time_point(now - start_time).count(); - int steps_per_second = elapsed > 0 ? iteration / elapsed : 0; + double steps_per_second = elapsed > 0 ? real_time(Decimal(iteration) * step_size) / elapsed : 0; - std::cout << "\r" << bar << " " << std::fixed << std::setprecision(2) - << percent << "% " << steps_per_second << " steps/s" << std::flush; + // 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 diff --git a/src/calc.hpp b/src/calc.hpp index ddcad40..054c322 100644 --- a/src/calc.hpp +++ b/src/calc.hpp @@ -12,7 +12,7 @@ std::vector> calculate_distances(const std::vector start_time, int length = 50); +void print_progress_bar(int iteration, int total, std::chrono::time_point start_time, int length, Decimal step_size); // Terminal plotting functions (if needed) #ifdef NCURSES_ENABLED diff --git a/src/main.cpp b/src/main.cpp index fb1edcc..7ea361b 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,41 +1,164 @@ -#include "simulator.hpp" -#include "units.hpp" #include +#include +#include +#include +#include +#include +#include "simulator.hpp" +#include "body.hpp" +#include "units.hpp" -int main() { +using json = nlohmann::json; +namespace po = boost::program_options; + +struct SimulationConfig { + std::string config_file; + std::string output_file; + int steps; + int steps_per_save; + Decimal step_size; // Changed to Decimal for normalized time + bool overwrite_output; + double simulation_time; // in seconds (real time) +}; + +// Convert time string to seconds +double parse_time(const std::string& time_str) { + std::string value_str = time_str; + std::string unit; + + // Extract the unit (last character) + if (!time_str.empty()) { + unit = time_str.back(); + value_str = time_str.substr(0, time_str.length() - 1); + } + + double value = std::stod(value_str); + + // Convert to seconds based on unit + switch (unit[0]) { + case 's': return value; // seconds + case 'm': return value * 60; // minutes + case 'h': return value * 3600; // hours + case 'd': return value * 86400; // days + default: throw std::runtime_error("Invalid time unit. Use s/m/h/d for seconds/minutes/hours/days"); + } +} + +SimulationConfig parse_command_line(int argc, char* argv[]) { + SimulationConfig config; + double temp_step_size; // Temporary variable for parsing + + po::options_description desc("Orbital Simulator Options"); + desc.add_options() + ("help,h", "Show help message") + ("config,c", po::value(&config.config_file)->required(), + "Path to planet configuration file (JSON)") + ("output,o", po::value(&config.output_file)->required(), + "Path to output file") + ("time,t", po::value()->required(), + "Simulation time with unit (e.g., 1h for 1 hour, 30m for 30 minutes, 2d for 2 days)") + ("step-size,s", po::value(&temp_step_size)->default_value(1.0), + "Simulation step size in seconds") + ("steps-per-save,p", po::value(&config.steps_per_save)->default_value(100), + "Number of steps between saves") + ("overwrite,w", po::bool_switch(&config.overwrite_output), + "Overwrite output file if it exists"); + + po::variables_map vm; try { - // Create bodies - std::vector bodies; + po::store(po::parse_command_line(argc, argv, desc), vm); - // 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"); + if (vm.count("help")) { + std::cout << desc << "\n"; + exit(0); + } - // 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"); + po::notify(vm); - // 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 + // Parse simulation time + config.simulation_time = parse_time(vm["time"].as()); - Simulator sim(bodies, step_size, steps_per_save, std::filesystem::path("output.txt")); + // Convert step size to Decimal and normalize + config.step_size = norm_time(Decimal(temp_step_size)); - // 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"; + // Calculate number of steps based on normalized time and step size + config.steps = static_cast(norm_time(config.simulation_time) / config.step_size); + } catch (const po::error& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << desc << "\n"; + exit(1); + } catch (const std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << desc << "\n"; + exit(1); + } + + return config; +} + +std::vector load_planets(const std::string& config_file) { + std::ifstream f(config_file); + if (!f.is_open()) { + throw std::runtime_error("Could not open config file: " + config_file); + } + + json j; + f >> j; + + std::vector bodies; + for (const auto& planet : j["planets"]) { + std::string name = planet["name"]; + double mass = planet["mass"]; + + std::vector pos = planet["position"]; + std::vector vel = planet["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 planets from config file + auto bodies = load_planets(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; diff --git a/src/simulator.cpp b/src/simulator.cpp index 67010f3..42ccb4d 100644 --- a/src/simulator.cpp +++ b/src/simulator.cpp @@ -29,7 +29,7 @@ Simulator::Simulator(const std::vector& bodies, } // Write initial header with masses - out = std::ofstream(output_file, std::ios::app); + out.open(output_file, std::ios::app); if (!out) { throw std::runtime_error("Failed to open output file: " + output_file.string()); } @@ -55,11 +55,9 @@ void Simulator::run(int steps) { 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); + if (i % steps_per_save == 0) { checkpoint(); + print_progress_bar(i + 1, steps, start_time, 50, step_size); } } } @@ -78,7 +76,7 @@ void Simulator::calculate_forces() { body.setAcceleration(Acceleration{Decimal(0), Decimal(0), Decimal(0)}); } - for (size_t i = 0; i < bodies.size(); i++) { + 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) { @@ -86,22 +84,22 @@ void Simulator::calculate_forces() { } // Calculate distance - Decimal dist = boost::multiprecision::sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); + 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 - Decimal force_magnitude = G * bodies[i].getMass() * bodies[j].getMass() / (dist * dist); + // F = G * m1 * m2 / r^2 BUT G = 1, and we'll multiply by the opposite mass later + // for the acceleration + Decimal force_magnitude = 1 / (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(); - + 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 / dist; // Negative because force is attractive - acc_j[k] = vec[k] * acc_magnitude_j / dist; // Positive because force is attractive + acc_i[k] = -vec[k] * acc_magnitude_i / dist; + acc_j[k] = vec[k] * acc_magnitude_j / dist; } // Add to current accelerations @@ -124,26 +122,9 @@ void Simulator::move_bodies() { } 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"; + out << body.toString() << "\n"; } + out << "\n"; + out.flush(); } \ No newline at end of file diff --git a/src/simulator.hpp b/src/simulator.hpp index 01ee2d0..58a8aeb 100644 --- a/src/simulator.hpp +++ b/src/simulator.hpp @@ -5,6 +5,7 @@ #include #include #include +#include class Simulator { public: diff --git a/src/units.hpp b/src/units.hpp index bb40a7f..2097234 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -14,24 +14,24 @@ using Acceleration = std::array; 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 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 = Decimal("734767309e14"); -const Decimal MOON_ORBITAL_VELOCITY = Decimal("1022"); // m/s relative to earth +const Decimal MOON_MASS = 734767309e14;//Decimal("734767309e14"); +const Decimal MOON_ORBITAL_VELOCITY = 1022;//Decimal("1022"); // m/s relative to earth -const Decimal SUN_MASS = Decimal("1989e27"); // kg -const Decimal SUN_RADIUS = Decimal("6957e5"); // meters +const Decimal SUN_MASS = 1989e27;//Decimal("1989e27"); // kg +const Decimal SUN_RADIUS = 6957e5;//Decimal("6957e5"); // meters -const Decimal PI = Decimal("3.14159265358979323846264338327950288419716939937510"); +const Decimal PI = 3.14159265358979323846264338327950288419716939937510; // Normalizing constants -const Decimal G = Decimal("6.67430e-11"); +const Decimal G = 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)); +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; }