Added config file for planets

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

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

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]
+        }
+    ]
+} 

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]
+        }
+    ]
+} 

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 {

src/calc.cpp

@@ -16,7 +16,7 @@ std::vector<std::vector<Decimal>> calculate_distances(const std::vector<Position
                 Decimal diff = positions[i][k] - positions[j][k];
                 sum += diff * diff;
             }
-            Decimal d = boost::multiprecision::sqrt(sum);
+            Decimal d = std::sqrt(sum);
             dists[i][j] = d;
             dists[j][i] = d;
         }
@@ -34,7 +34,7 @@ std::string format_sig_figs(Decimal value, int sig_figs) {
     return ss.str();
 }
 
-void print_progress_bar(int iteration, int total, std::chrono::time_point<std::chrono::steady_clock> start_time, int length) {
+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;
     
@@ -47,10 +47,26 @@ void print_progress_bar(int iteration, int total, std::chrono::time_point<std::c
     
     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;
+    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

src/calc.hpp

@@ -12,7 +12,7 @@ std::vector<std::vector<Decimal>> calculate_distances(const std::vector<Position
 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);
+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

src/main.cpp

@@ -1,41 +1,164 @@
-#include "simulator.hpp"
-#include "units.hpp"
 #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"
 
-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<std::string>(&config.config_file)->required(), 
+         "Path to planet 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 {
-        // Create bodies
+        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_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<Body> bodies;
+    for (const auto& planet : j["planets"]) {
+        std::string name = planet["name"];
+        double mass = planet["mass"];
         
-        // 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");
+        std::vector<double> pos = planet["position"];
+        std::vector<double> vel = planet["velocity"];
         
-        // 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");
+        if (pos.size() != 3 || vel.size() != 3) {
+            throw std::runtime_error("Position and velocity must be 3D vectors");
+        }
         
-        // 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
+        // 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);
         
-        Simulator sim(bodies, step_size, steps_per_save, std::filesystem::path("output.txt"));
+        bodies.emplace_back(position, velocity, normalized_mass, name);
+        std::cout << "Loaded " << name << " with mass " << mass << " kg\n";
+    }
     
-        // 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 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;

src/simulator.cpp

@@ -29,7 +29,7 @@ Simulator::Simulator(const std::vector<Body>& 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);
         }
     }
 }
@@ -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 << body.toString() << "\n";
     }
     out << "\n";
-    }
+    out.flush();
 } 

src/simulator.hpp

@@ -5,6 +5,7 @@
 #include <string>
 #include <filesystem>
 #include <memory>
+#include <fstream>
 
 class Simulator {
 public:

src/units.hpp

@@ -14,24 +14,24 @@ 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 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; }