Added more planets, plotting, json AND toml support

PLOTTING.md

@@ -0,0 +1,88 @@
+# Trajectory Plotting
+
+This Python script can read and visualize the binary trajectory files generated by the orbital simulator.
+
+## Installation
+
+First, install the required Python dependencies:
+
+```bash
+pip install -r requirements.txt
+```
+
+## Usage
+
+### Basic Usage
+
+```bash
+# Plot trajectories from a simulation output file
+python plot_trajectories.py trajectory_output.bin
+```
+
+### Options
+
+- `--2d-only`: Only show 2D plot (X-Y plane)  
+- `--3d-only`: Only show 3D plot
+- `--energy`: Also show energy conservation plot
+- `--animate`: Show animated trajectories with real-time simulation
+- `--save-animation <filename>`: Save animation as MP4 file (requires ffmpeg)
+- `--static`: Show static plots (default behavior)
+- `--interval <ms>`: Animation frame interval in milliseconds (default: auto-scaled)
+- `--target-duration <seconds>`: Target animation duration (default: 60)
+- `--center <body_name>`: Center animation/plot on specified body (e.g., "Sun", "Earth")
+- `--list-bodies`: List available bodies in trajectory file and exit
+
+### Examples
+
+```bash
+# Run a simulation and plot the results
+cargo run --bin simulator -- --config config/planets.json --time 365d --step-size 3600 --output-file solar_system.bin
+python plot_trajectories.py solar_system.bin
+
+# List available bodies in trajectory file
+python plot_trajectories.py solar_system.bin --list-bodies
+
+# Show animated trajectories centered on the Sun
+python plot_trajectories.py solar_system.bin --animate --center Sun
+
+# Show Earth-centered view (great for seeing Moon's orbit)
+python plot_trajectories.py solar_system.bin --animate --center Earth --2d-only
+
+# Show only 2D animation with custom speed
+python plot_trajectories.py solar_system.bin --animate --2d-only --interval 100
+
+# Save animation as MP4 file (Sun-centered)
+python plot_trajectories.py solar_system.bin --animate --center Sun --save-animation solar_system_animation
+
+# Show 3D trajectories with energy plot (static, Earth-centered)
+python plot_trajectories.py solar_system.bin --3d-only --energy --static --center Earth
+```
+
+## Output
+
+The script can display:
+
+### Static Plots
+- **2D Plot**: Trajectories in the X-Y plane with starting positions (circles) and ending positions (squares)
+- **3D Plot**: Full 3D orbital trajectories  
+- **Energy Plot** (optional): Shows kinetic and total energy over time to verify energy conservation
+
+### Animated Plots
+- **2D Animation**: Real-time orbital motion in the X-Y plane with time display
+- **3D Animation**: Full 3D animated orbital trajectories
+- **Time Display**: Shows current simulation time in seconds/hours/days
+- **Trails**: Recent positions shown as fading trails behind each body
+- **MP4 Export**: Save animations as video files (requires ffmpeg)
+
+Each celestial body is plotted in a different color with its name in the legend.
+
+## File Format
+
+The binary trajectory file contains serialized snapshots with:
+- Timestamp (f64)
+- Array of bodies, each containing:
+  - Name (String)
+  - Mass (f64) 
+  - Position (3x f64)
+  - Velocity (3x f64)
+  - Acceleration (3x f64)

config/jwst.json

@@ -1,22 +0,0 @@
-{
-    "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]
-        }
-    ]
-} 

config/planets.json

@@ -1,34 +0,0 @@
-{
-    "bodies": [
-        {
-            "name": "Mercury",
-            "mass": 3.30104e23,
-            "position": [4.6000e10, 0, 0],
-            "velocity": [0, 58970, 0]
-        },
-        {
-            "name": "Venus",
-            "mass": 4.867e24,
-            "position": [1.08941e11, 0, 0],
-            "velocity": [0, 34780, 0]
-        },
-        {
-            "name": "Earth",
-            "mass": 5.972e24,
-            "position": [1.47095e11, 0, 0],
-            "velocity": [0, 29290, 0]
-        },
-        {
-            "name": "Moon",
-            "mass": 7.34767309e22,
-            "position": [149982270700, 0, 0],
-            "velocity": [0, 30822, 0]
-        },
-        {
-            "name": "Sun",
-            "mass": 1.989e30,
-            "position": [0, 0, 0],
-            "velocity": [0, 0, 0]
-        }
-    ]
-} 

config/planets.toml

@@ -1,9 +1,16 @@
 # Simulation body configuration
+
+[[bodies]]
+name = "Sun"
+mass = 1.989e30
+position = [0.0, 0.0, 0.0]
+velocity = [0.0, 0.0, 0.0]
+
 [[bodies]]
 name = "Mercury"
-mass = 3.30104e23         # kg
+mass = 3.30104e23
-position = [46000000000.0, 0.0, 0.0]  # meters
+position = [4.6000e10, 0.0, 0.0]        # 0.307 AU
-velocity = [0.0, 58970.0, 0.0]        # m/s
+velocity = [0.0, 58970.0, 0.0]          # m/s
 
 [[bodies]]
 name = "Venus"
@@ -15,16 +22,80 @@ velocity = [0.0, 34780.0, 0.0]
 name = "Earth"
 mass = 5.972e24
 position = [147095000000.0, 0.0, 0.0]
-velocity = [0.0, 29290.0, 0.0]
+velocity = [0.0, 30290.0, 0.0]
 
 [[bodies]]
 name = "Moon"
 mass = 7.34767309e22
-position = [149982270700.0, 0.0, 0.0]
+position = [147458300000, 0.0, 0.0]       # Earth + 384,400 km
-velocity = [0.0, 30822.0, 0.0]
+velocity = [0.0, 31372.0, 0.0]          # Earth velocity + moon orbital velocity
 
 [[bodies]]
-name = "Sun"
+name = "Mars"
-mass = 1.989e30
+mass = 6.4171e23
-position = [0.0, 0.0, 0.0]
+position = [2.279e11, 0.0, 0.0]         # 1.524 AU
-velocity = [0.0, 0.0, 0.0]
+velocity = [0.0, 24007.0, 0.0]          # m/s
+
+[[bodies]]
+name = "Jupiter"
+mass = 1.8982e27
+position = [7.786e11, 0.0, 0.0]         # 5.204 AU
+velocity = [0.0, 13060.0, 0.0]          # m/s
+
+[[bodies]]
+name = "Saturn"
+mass = 5.6834e26
+position = [1.432e12, 0.0, 0.0]         # 9.573 AU
+velocity = [0.0, 9640.0, 0.0]           # m/s
+
+[[bodies]]
+name = "Uranus"
+mass = 8.6810e25
+position = [2.867e12, 0.0, 0.0]         # 19.165 AU
+velocity = [0.0, 6810.0, 0.0]           # m/s
+
+[[bodies]]
+name = "Neptune"
+mass = 1.02413e26
+position = [4.515e12, 0.0, 0.0]         # 30.178 AU
+velocity = [0.0, 5430.0, 0.0]           # m/s
+
+# Dwarf planets and interesting objects
+
+[[bodies]]
+name = "Pluto"
+mass = 1.303e22
+position = [5.906e12, 0.0, 0.0]         # 39.482 AU (average)
+velocity = [0.0, 4670.0, 0.0]           # m/s
+
+[[bodies]]
+name = "Ceres"
+mass = 9.393e20
+position = [4.14e11, 0.0, 0.0]          # 2.766 AU (asteroid belt)
+velocity = [0.0, 17880.0, 0.0]          # m/s
+
+# Some major moons for more interesting dynamics
+
+[[bodies]]
+name = "Io"
+mass = 8.932e22
+position = [7.790e11, 0.0, 0.0]         # Jupiter + 421,700 km
+velocity = [0.0, 30350.0, 0.0]          # Jupiter velocity + Io orbital velocity
+
+[[bodies]]
+name = "Europa"
+mass = 4.800e22
+position = [7.793e11, 0.0, 0.0]         # Jupiter + 671,034 km
+velocity = [0.0, 26890.0, 0.0]          # Jupiter velocity + Europa orbital velocity
+
+[[bodies]]
+name = "Ganymede"
+mass = 1.482e23
+position = [7.796e11, 0.0, 0.0]         # Jupiter + 1,070,412 km
+velocity = [0.0, 23250.0, 0.0]          # Jupiter velocity + Ganymede orbital velocity
+
+[[bodies]]
+name = "Titan"
+mass = 1.345e23
+position = [1.433e12, 0.0, 0.0]         # Saturn + 1,221,830 km
+velocity = [0.0, 15100.0, 0.0]          # Saturn velocity + Titan orbital velocity

config/planets_fake.json

@@ -1,16 +0,0 @@
-{
-    "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]
-        }
-    ]
-} 

inspect_trajectories.py

@@ -0,0 +1,347 @@
+#!/usr/bin/env python3
+"""
+Inspect orbital trajectory data from binary files in a tabular format.
+
+This script reads the binary trajectory files generated by the orbital simulator
+and displays the data in a nicely formatted table for inspection and debugging.
+
+Usage:
+    python inspect_trajectories.py <trajectory_file.bin>
+"""
+
+import sys
+import struct
+import argparse
+import numpy as np
+from collections import defaultdict
+
+class BinaryReader:
+    def __init__(self, data):
+        self.data = data
+        self.pos = 0
+    
+    def read_u64(self):
+        result = struct.unpack('<Q', self.data[self.pos:self.pos+8])[0]
+        self.pos += 8
+        return result
+    
+    def read_f64(self):
+        result = struct.unpack('<d', self.data[self.pos:self.pos+8])[0]
+        self.pos += 8
+        return result
+    
+    def read_string(self):
+        # Read length (u64) then string bytes
+        length = self.read_u64()
+        result = self.data[self.pos:self.pos+length].decode('utf-8')
+        self.pos += length
+        return result
+    
+    def read_vec3(self):
+        # Read 3 f64 values for position/velocity/acceleration
+        x = self.read_f64()
+        y = self.read_f64()
+        z = self.read_f64()
+        return np.array([x, y, z])
+
+def read_trajectory_file(filename):
+    """Read the binary trajectory file and return parsed data."""
+    with open(filename, 'rb') as f:
+        data = f.read()
+    
+    reader = BinaryReader(data)
+    snapshots = []
+    
+    try:
+        while reader.pos < len(data):
+            # Read snapshot
+            time = reader.read_f64()
+            
+            # Read number of bodies (u64)
+            num_bodies = reader.read_u64()
+            
+            bodies = []
+            for _ in range(num_bodies):
+                # Read Body struct
+                name = reader.read_string()
+                mass = reader.read_f64()
+                position = reader.read_vec3()
+                velocity = reader.read_vec3()
+                acceleration = reader.read_vec3()
+                
+                bodies.append({
+                    'name': name,
+                    'mass': mass,
+                    'position': position,
+                    'velocity': velocity,
+                    'acceleration': acceleration
+                })
+            
+            snapshots.append({
+                'time': time,
+                'bodies': bodies
+            })
+            
+    except struct.error:
+        # End of file or corrupted data
+        pass
+    
+    return snapshots
+
+def format_scientific(value, precision=3):
+    """Format number in scientific notation with specified precision."""
+    if abs(value) < 1e-6 or abs(value) >= 1e6:
+        return f"{value:.{precision}e}"
+    else:
+        return f"{value:.{precision}f}"
+
+def format_vector(vec, precision=3):
+    """Format a 3D vector in a compact form."""
+    x, y, z = vec
+    return f"({format_scientific(x, precision)}, {format_scientific(y, precision)}, {format_scientific(z, precision)})"
+
+def format_time(seconds):
+    """Format time in a human-readable format."""
+    if seconds == 0:
+        return "0.0s"
+    elif seconds < 60:
+        return f"{seconds:.1f}s"
+    elif seconds < 3600:
+        return f"{seconds/60:.1f}m"
+    elif seconds < 86400:
+        return f"{seconds/3600:.1f}h"
+    else:
+        return f"{seconds/86400:.1f}d"
+
+def print_summary_table(snapshots, max_rows=15):
+    """Print a summary table of all snapshots."""
+    print("📊 TRAJECTORY SUMMARY")
+    print("=" * 100)
+    
+    # Header
+    print(f"{'Step':<6} {'Time':<12} {'Bodies':<8} {'Sample Position (first body)':<35} {'Sample Velocity (first body)':<35}")
+    print("-" * 100)
+    
+    # Determine which snapshots to show
+    total_snapshots = len(snapshots)
+    if total_snapshots <= max_rows:
+        indices = list(range(total_snapshots))
+    else:
+        # Show first few, middle few, and last few
+        start_count = max_rows // 3
+        end_count = max_rows // 3
+        middle_count = max_rows - start_count - end_count
+        
+        indices = []
+        indices.extend(range(start_count))
+        
+        if middle_count > 0:
+            middle_start = total_snapshots // 2 - middle_count // 2
+            indices.extend(range(middle_start, middle_start + middle_count))
+        
+        indices.extend(range(total_snapshots - end_count, total_snapshots))
+        indices = sorted(set(indices))  # Remove duplicates and sort
+    
+    prev_index = -1
+    for i, idx in enumerate(indices):
+        if idx > prev_index + 1 and prev_index >= 0:
+            print("    ...")
+        
+        snapshot = snapshots[idx]
+        time_str = format_time(snapshot['time'])
+        body_count = len(snapshot['bodies'])
+        
+        if snapshot['bodies']:
+            first_body = snapshot['bodies'][0]
+            pos_str = format_vector(first_body['position'], 2)
+            vel_str = format_vector(first_body['velocity'], 2)
+        else:
+            pos_str = "N/A"
+            vel_str = "N/A"
+        
+        print(f"{idx:<6} {time_str:<12} {body_count:<8} {pos_str:<35} {vel_str:<35}")
+        prev_index = idx
+    
+    print("-" * 100)
+    print(f"Total snapshots: {total_snapshots}")
+
+def print_detailed_table(snapshots, body_name=None, max_rows=15):
+    """Print detailed information for a specific body."""
+    if not snapshots:
+        print("No data to display!")
+        return
+    
+    # Get all body names
+    all_bodies = set()
+    for snapshot in snapshots:
+        for body in snapshot['bodies']:
+            all_bodies.add(body['name'])
+    
+    if body_name and body_name not in all_bodies:
+        print(f"❌ Body '{body_name}' not found. Available bodies: {', '.join(sorted(all_bodies))}")
+        return
+    
+    # If no body specified, show the first one
+    if not body_name:
+        body_name = sorted(all_bodies)[0]
+    
+    print(f"🌍 DETAILED VIEW: {body_name}")
+    print("=" * 120)
+    
+    # Header
+    print(f"{'Step':<6} {'Time':<12} {'Position (x, y, z)':<40} {'Velocity (x, y, z)':<40} {'|v|':<12} {'|a|':<12}")
+    print("-" * 120)
+    
+    # Determine which snapshots to show
+    total_snapshots = len(snapshots)
+    if total_snapshots <= max_rows:
+        indices = list(range(total_snapshots))
+    else:
+        # Show first few, middle few, and last few
+        start_count = max_rows // 3
+        end_count = max_rows // 3
+        middle_count = max_rows - start_count - end_count
+        
+        indices = []
+        indices.extend(range(start_count))
+        
+        if middle_count > 0:
+            middle_start = total_snapshots // 2 - middle_count // 2
+            indices.extend(range(middle_start, middle_start + middle_count))
+        
+        indices.extend(range(total_snapshots - end_count, total_snapshots))
+        indices = sorted(set(indices))
+    
+    prev_index = -1
+    for idx in indices:
+        if idx > prev_index + 1 and prev_index >= 0:
+            print("    ...")
+        
+        snapshot = snapshots[idx]
+        time_str = format_time(snapshot['time'])
+        
+        # Find the body in this snapshot
+        body_data = None
+        for body in snapshot['bodies']:
+            if body['name'] == body_name:
+                body_data = body
+                break
+        
+        if body_data:
+            pos_str = format_vector(body_data['position'], 3)
+            vel_str = format_vector(body_data['velocity'], 3)
+            vel_mag = np.linalg.norm(body_data['velocity'])
+            acc_mag = np.linalg.norm(body_data['acceleration'])
+            vel_mag_str = format_scientific(vel_mag, 2)
+            acc_mag_str = format_scientific(acc_mag, 2)
+            
+            print(f"{idx:<6} {time_str:<12} {pos_str:<40} {vel_str:<40} {vel_mag_str:<12} {acc_mag_str:<12}")
+        else:
+            print(f"{idx:<6} {time_str:<12} {'BODY NOT FOUND':<40}")
+        
+        prev_index = idx
+    
+    print("-" * 120)
+
+def print_statistics(snapshots):
+    """Print statistical information about the trajectory."""
+    if not snapshots:
+        return
+    
+    print("\n📈 TRAJECTORY STATISTICS")
+    print("=" * 80)
+    
+    # Time statistics
+    times = [s['time'] for s in snapshots]
+    time_start, time_end = times[0], times[-1]
+    duration = time_end - time_start
+    
+    print(f"Time range: {format_time(time_start)} to {format_time(time_end)}")
+    print(f"Total duration: {format_time(duration)}")
+    print(f"Number of snapshots: {len(snapshots)}")
+    
+    if len(times) > 1:
+        time_steps = [times[i+1] - times[i] for i in range(len(times)-1)]
+        avg_step = np.mean(time_steps)
+        print(f"Average time step: {format_time(avg_step)}")
+    
+    # Body statistics
+    body_names = set()
+    for snapshot in snapshots:
+        for body in snapshot['bodies']:
+            body_names.add(body['name'])
+    
+    print(f"Bodies tracked: {', '.join(sorted(body_names))}")
+    
+    # Position and velocity ranges for each body
+    for body_name in sorted(body_names):
+        positions = []
+        velocities = []
+        
+        for snapshot in snapshots:
+            for body in snapshot['bodies']:
+                if body['name'] == body_name:
+                    positions.append(body['position'])
+                    velocities.append(body['velocity'])
+        
+        if positions:
+            positions = np.array(positions)
+            velocities = np.array(velocities)
+            
+            pos_min = np.min(positions, axis=0)
+            pos_max = np.max(positions, axis=0)
+            vel_min = np.min(np.linalg.norm(velocities, axis=1))
+            vel_max = np.max(np.linalg.norm(velocities, axis=1))
+            
+            print(f"\n{body_name}:")
+            print(f"  Position range: X[{format_scientific(pos_min[0])}, {format_scientific(pos_max[0])}]")
+            print(f"                  Y[{format_scientific(pos_min[1])}, {format_scientific(pos_max[1])}]")
+            print(f"                  Z[{format_scientific(pos_min[2])}, {format_scientific(pos_max[2])}]")
+            print(f"  Speed range: {format_scientific(vel_min)} to {format_scientific(vel_max)} m/s")
+
+def main():
+    parser = argparse.ArgumentParser(description='Inspect orbital trajectory data in tabular format')
+    parser.add_argument('trajectory_file', help='Binary trajectory file to inspect')
+    parser.add_argument('--rows', '-r', type=int, default=15, help='Maximum number of rows to display (default: 15)')
+    parser.add_argument('--body', '-b', type=str, help='Show detailed view for specific body')
+    parser.add_argument('--summary', '-s', action='store_true', help='Show summary of all snapshots')
+    parser.add_argument('--stats', action='store_true', help='Show trajectory statistics')
+    parser.add_argument('--all', '-a', action='store_true', help='Show all available information')
+    
+    args = parser.parse_args()
+    
+    print(f"🔍 Inspecting trajectory file: {args.trajectory_file}")
+    print()
+    
+    try:
+        snapshots = read_trajectory_file(args.trajectory_file)
+        
+        if not snapshots:
+            print("❌ No data found in file!")
+            return
+        
+        # Determine what to show
+        show_summary = args.summary or args.all or (not args.body and not args.stats)
+        show_detailed = args.body or args.all
+        show_stats = args.stats or args.all
+        
+        if show_summary:
+            print_summary_table(snapshots, args.rows)
+            print()
+        
+        if show_detailed:
+            print_detailed_table(snapshots, args.body, args.rows)
+            print()
+        
+        if show_stats:
+            print_statistics(snapshots)
+        
+    except FileNotFoundError:
+        print(f"❌ Error: File '{args.trajectory_file}' not found!")
+    except Exception as e:
+        print(f"❌ Error reading file: {e}")
+        import traceback
+        traceback.print_exc()
+
+if __name__ == "__main__":
+    main()

plot_trajectories.py

@@ -0,0 +1,570 @@
+#!/usr/bin/env python3
+"""
+Plot orbital trajectories from binary output file generated by the orbital simulator.
+
+The binary file contains snapshots serialized with bincode, where each snapshot has:
+- time: f64
+- bodies: array of Body structs
+  - name: String
+  - mass: f64
+  - position: [f64; 3] (x, y, z)
+  - velocity: [f64; 3] (vx, vy, vz)
+  - acceleration: [f64; 3] (ax, ay, az)
+
+Usage:
+    python plot_trajectories.py <trajectory_file.bin>
+"""
+
+import sys
+import struct
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.animation as animation
+import numpy as np
+from collections import defaultdict
+import argparse
+
+class BinaryReader:
+    def __init__(self, data):
+        self.data = data
+        self.pos = 0
+    
+    def read_u64(self):
+        result = struct.unpack('<Q', self.data[self.pos:self.pos+8])[0]
+        self.pos += 8
+        return result
+    
+    def read_f64(self):
+        result = struct.unpack('<d', self.data[self.pos:self.pos+8])[0]
+        self.pos += 8
+        return result
+    
+    def read_string(self):
+        # Read length (u64) then string bytes
+        length = self.read_u64()
+        result = self.data[self.pos:self.pos+length].decode('utf-8')
+        self.pos += length
+        return result
+    
+    def read_vec3(self):
+        # Read 3 f64 values for position/velocity/acceleration
+        x = self.read_f64()
+        y = self.read_f64()
+        z = self.read_f64()
+        return np.array([x, y, z])
+
+def read_trajectory_file(filename):
+    """Read the binary trajectory file and return parsed data."""
+    with open(filename, 'rb') as f:
+        data = f.read()
+    
+    reader = BinaryReader(data)
+    snapshots = []
+    
+    try:
+        while reader.pos < len(data):
+            # Read snapshot
+            time = reader.read_f64()
+            
+            # Read number of bodies (u64)
+            num_bodies = reader.read_u64()
+            
+            bodies = []
+            for _ in range(num_bodies):
+                # Read Body struct
+                name = reader.read_string()
+                mass = reader.read_f64()
+                position = reader.read_vec3()
+                velocity = reader.read_vec3()
+                acceleration = reader.read_vec3()
+                
+                bodies.append({
+                    'name': name,
+                    'mass': mass,
+                    'position': position,
+                    'velocity': velocity,
+                    'acceleration': acceleration
+                })
+            
+            snapshots.append({
+                'time': time,
+                'bodies': bodies
+            })
+            
+    except struct.error:
+        # End of file or corrupted data
+        pass
+    
+    return snapshots
+
+def organize_trajectories(snapshots):
+    """Organize snapshots into trajectories by body name."""
+    trajectories = defaultdict(list)
+    times = []
+    
+    for snapshot in snapshots:
+        times.append(snapshot['time'])
+        for body in snapshot['bodies']:
+            trajectories[body['name']].append(body['position'])
+    
+    # Convert lists to numpy arrays
+    for name in trajectories:
+        trajectories[name] = np.array(trajectories[name])
+    
+    return dict(trajectories), np.array(times)
+
+def plot_trajectories_2d(trajectories, times, center_body=None):
+    """Plot 2D trajectories (X-Y plane)."""
+    plt.figure(figsize=(12, 10))
+    
+    colors = plt.cm.tab10(np.linspace(0, 1, len(trajectories)))
+    
+    for i, (name, positions) in enumerate(trajectories.items()):
+        x = positions[:, 0]
+        y = positions[:, 1]
+        
+        plt.plot(x, y, color=colors[i], alpha=0.7, linewidth=1.5, label=name)
+        
+        # Mark starting position
+        plt.plot(x[0], y[0], 'o', color=colors[i], markersize=8, alpha=0.8)
+        
+        # Mark current position
+        plt.plot(x[-1], y[-1], 's', color=colors[i], markersize=6, alpha=0.8)
+    
+    plt.xlabel('X Position (m)')
+    plt.ylabel('Y Position (m)')
+    
+    title = 'Orbital Trajectories (X-Y Plane)'
+    if center_body:
+        title += f' - Centered on {center_body}'
+    plt.title(title)
+    
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.axis('equal')
+    plt.tight_layout()
+
+def plot_trajectories_3d(trajectories, times, center_body=None):
+    """Plot 3D trajectories."""
+    fig = plt.figure(figsize=(12, 10))
+    ax = fig.add_subplot(111, projection='3d')
+    
+    colors = plt.cm.tab10(np.linspace(0, 1, len(trajectories)))
+    
+    for i, (name, positions) in enumerate(trajectories.items()):
+        x = positions[:, 0]
+        y = positions[:, 1]
+        z = positions[:, 2]
+        
+        ax.plot(x, y, z, color=colors[i], alpha=0.7, linewidth=1.5, label=name)
+        
+        # Mark starting position
+        ax.scatter(x[0], y[0], z[0], color=colors[i], s=100, alpha=0.8, marker='o')
+        
+        # Mark current position
+        ax.scatter(x[-1], y[-1], z[-1], color=colors[i], s=60, alpha=0.8, marker='s')
+    
+    ax.set_xlabel('X Position (m)')
+    ax.set_ylabel('Y Position (m)')
+    ax.set_zlabel('Z Position (m)')
+    
+    title = 'Orbital Trajectories (3D)'
+    if center_body:
+        title += f' - Centered on {center_body}'
+    ax.set_title(title)
+    
+    ax.legend()
+    
+    # Make axes equal
+    max_range = 0
+    for positions in trajectories.values():
+        range_val = np.max(np.abs(positions))
+        max_range = max(max_range, range_val)
+    
+    ax.set_xlim(-max_range, max_range)
+    ax.set_ylim(-max_range, max_range)
+    ax.set_zlim(-max_range, max_range)
+
+def plot_energy_over_time(snapshots, times):
+    """Plot energy evolution over time."""
+    plt.figure(figsize=(12, 6))
+    
+    total_energies = []
+    kinetic_energies = []
+    
+    for snapshot in snapshots:
+        ke = 0
+        pe = 0
+        bodies = snapshot['bodies']
+        
+        # Calculate kinetic energy
+        for body in bodies:
+            v_squared = np.sum(body['velocity']**2)
+            ke += 0.5 * body['mass'] * v_squared
+        
+        # Calculate potential energy (simplified, assuming G=1 in normalized units)
+        G = 6.67430e-11  # You might need to adjust this based on your normalization
+        for i in range(len(bodies)):
+            for j in range(i+1, len(bodies)):
+                r = np.linalg.norm(bodies[i]['position'] - bodies[j]['position'])
+                pe -= G * bodies[i]['mass'] * bodies[j]['mass'] / r
+        
+        kinetic_energies.append(ke)
+        total_energies.append(ke + pe)
+    
+    plt.plot(times, kinetic_energies, label='Kinetic Energy', alpha=0.8)
+    plt.plot(times, total_energies, label='Total Energy', alpha=0.8)
+    
+    plt.xlabel('Time (s)')
+    plt.ylabel('Energy (J)')
+    plt.title('Energy Conservation Over Time')
+    plt.legend()
+    plt.grid(True, alpha=0.3)
+    plt.tight_layout()
+
+def plot_animated_2d(trajectories, times, interval=50, center_body=None):
+    """Create animated 2D trajectory plot."""
+    fig, ax = plt.subplots(figsize=(12, 10))
+    
+    colors = plt.cm.tab10(np.linspace(0, 1, len(trajectories)))
+    body_names = list(trajectories.keys())
+    
+    # Set up the plot
+    ax.set_xlabel('X Position (m)')
+    ax.set_ylabel('Y Position (m)')
+    
+    title = 'Animated Orbital Trajectories (X-Y Plane)'
+    if center_body:
+        title += f' - Centered on {center_body}'
+    ax.set_title(title)
+    ax.grid(True, alpha=0.3)
+    
+    # Calculate plot limits
+    all_positions = np.concatenate(list(trajectories.values()))
+    margin = 0.1
+    x_range = np.max(all_positions[:, 0]) - np.min(all_positions[:, 0])
+    y_range = np.max(all_positions[:, 1]) - np.min(all_positions[:, 1])
+    x_center = (np.max(all_positions[:, 0]) + np.min(all_positions[:, 0])) / 2
+    y_center = (np.max(all_positions[:, 1]) + np.min(all_positions[:, 1])) / 2
+    
+    max_range = max(x_range, y_range) * (1 + margin)
+    ax.set_xlim(x_center - max_range/2, x_center + max_range/2)
+    ax.set_ylim(y_center - max_range/2, y_center + max_range/2)
+    ax.set_aspect('equal')
+    
+    # Initialize plot elements
+    trajectory_lines = []
+    body_points = []
+    body_trails = []
+    
+    for i, name in enumerate(body_names):
+        # Trajectory line (will grow over time)
+        line, = ax.plot([], [], color=colors[i], alpha=0.7, linewidth=1.5, label=name)
+        trajectory_lines.append(line)
+        
+        # Current body position
+        point, = ax.plot([], [], 'o', color=colors[i], markersize=8, alpha=0.9)
+        body_points.append(point)
+        
+        # Trail of recent positions
+        trail, = ax.plot([], [], 'o', color=colors[i], markersize=3, alpha=0.3)
+        body_trails.append(trail)
+    
+    # Time display
+    time_text = ax.text(0.02, 0.98, '', transform=ax.transAxes, fontsize=12,
+                       verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
+    
+    ax.legend(loc='upper right')
+    
+    def animate(frame):
+        current_time = times[frame]
+        
+        # Format time for display
+        if current_time >= 86400:  # More than a day
+            time_str = f"Time: {current_time/86400:.1f} days"
+        elif current_time >= 3600:  # More than an hour
+            time_str = f"Time: {current_time/3600:.1f} hours"
+        else:
+            time_str = f"Time: {current_time:.1f} seconds"
+        
+        time_text.set_text(time_str)
+        
+        # Update each body
+        for i, name in enumerate(body_names):
+            positions = trajectories[name]
+            
+            # Update trajectory line (show path up to current time)
+            x_data = positions[:frame+1, 0]
+            y_data = positions[:frame+1, 1]
+            trajectory_lines[i].set_data(x_data, y_data)
+            
+            # Update current position
+            if frame < len(positions):
+                current_pos = positions[frame]
+                body_points[i].set_data([current_pos[0]], [current_pos[1]])
+                
+                # Update trail (last 20 positions)
+                trail_start = max(0, frame - 20)
+                trail_x = positions[trail_start:frame, 0]
+                trail_y = positions[trail_start:frame, 1]
+                body_trails[i].set_data(trail_x, trail_y)
+        
+        return trajectory_lines + body_points + body_trails + [time_text]
+    
+    num_frames = len(times)
+    anim = animation.FuncAnimation(fig, animate, frames=num_frames, 
+                                 interval=interval, blit=True, repeat=True)
+    
+    return fig, anim
+
+def plot_animated_3d(trajectories, times, interval=50, center_body=None):
+    """Create animated 3D trajectory plot."""
+    fig = plt.figure(figsize=(12, 10))
+    ax = fig.add_subplot(111, projection='3d')
+    
+    colors = plt.cm.tab10(np.linspace(0, 1, len(trajectories)))
+    body_names = list(trajectories.keys())
+    
+    # Set up the plot
+    ax.set_xlabel('X Position (m)')
+    ax.set_ylabel('Y Position (m)')
+    ax.set_zlabel('Z Position (m)')
+    
+    title = 'Animated Orbital Trajectories (3D)'
+    if center_body:
+        title += f' - Centered on {center_body}'
+    ax.set_title(title)
+    
+    # Calculate plot limits
+    all_positions = np.concatenate(list(trajectories.values()))
+    max_range = np.max(np.abs(all_positions)) * 1.1
+    ax.set_xlim(-max_range, max_range)
+    ax.set_ylim(-max_range, max_range)
+    ax.set_zlim(-max_range, max_range)
+    
+    # Initialize plot elements
+    trajectory_lines = []
+    body_points = []
+    
+    for i, name in enumerate(body_names):
+        # Trajectory line
+        line, = ax.plot([], [], [], color=colors[i], alpha=0.7, linewidth=1.5, label=name)
+        trajectory_lines.append(line)
+        
+        # Current body position
+        point = ax.scatter([], [], [], color=colors[i], s=100, alpha=0.9)
+        body_points.append(point)
+    
+    # Time display
+    time_text = ax.text2D(0.02, 0.98, '', transform=ax.transAxes, fontsize=12,
+                         verticalalignment='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
+    
+    ax.legend(loc='upper right')
+    
+    def animate(frame):
+        current_time = times[frame]
+        
+        # Format time for display
+        if current_time >= 86400:  # More than a day
+            time_str = f"Time: {current_time/86400:.1f} days"
+        elif current_time >= 3600:  # More than an hour
+            time_str = f"Time: {current_time/3600:.1f} hours"
+        else:
+            time_str = f"Time: {current_time:.1f} seconds"
+        
+        time_text.set_text(time_str)
+        
+        # Update each body
+        for i, name in enumerate(body_names):
+            positions = trajectories[name]
+            
+            # Update trajectory line
+            x_data = positions[:frame+1, 0]
+            y_data = positions[:frame+1, 1]
+            z_data = positions[:frame+1, 2]
+            trajectory_lines[i].set_data(x_data, y_data)
+            trajectory_lines[i].set_3d_properties(z_data)
+            
+            # Update current position
+            if frame < len(positions):
+                current_pos = positions[frame]
+                # Remove old scatter point and create new one
+                body_points[i].remove()
+                body_points[i] = ax.scatter([current_pos[0]], [current_pos[1]], [current_pos[2]], 
+                                         color=colors[i], s=100, alpha=0.9)
+        
+        return trajectory_lines + body_points + [time_text]
+    
+    num_frames = len(times)
+    anim = animation.FuncAnimation(fig, animate, frames=num_frames, 
+                                 interval=interval, blit=False, repeat=True)
+    
+    return fig, anim
+
+def calculate_animation_params(num_frames, target_duration_sec=60, manual_interval=None):
+    """Calculate animation parameters for optimal viewing experience."""
+    if manual_interval is not None:
+        # User specified manual interval
+        interval_ms = manual_interval
+        total_duration_sec = num_frames * interval_ms / 1000.0
+        time_scale_factor = target_duration_sec / total_duration_sec
+        return interval_ms, total_duration_sec, time_scale_factor, True
+    
+    # Auto-calculate interval for target duration
+    target_duration_ms = target_duration_sec * 1000
+    optimal_interval = max(10, target_duration_ms // num_frames)  # Minimum 10ms for smooth animation
+    actual_duration_sec = num_frames * optimal_interval / 1000.0
+    time_scale_factor = target_duration_sec / actual_duration_sec
+    
+    return optimal_interval, actual_duration_sec, time_scale_factor, False
+
+def center_trajectories_on_body(trajectories, center_body_name):
+    """Center all trajectories relative to the specified body."""
+    if center_body_name not in trajectories:
+        available_bodies = list(trajectories.keys())
+        raise ValueError(f"Body '{center_body_name}' not found. Available bodies: {available_bodies}")
+    
+    center_trajectory = trajectories[center_body_name]
+    centered_trajectories = {}
+    
+    for body_name, trajectory in trajectories.items():
+        # Subtract the center body's position from each body's trajectory
+        centered_trajectories[body_name] = trajectory - center_trajectory
+    
+    return centered_trajectories
+
+def main():
+    parser = argparse.ArgumentParser(description='Plot orbital trajectories from binary file')
+    parser.add_argument('trajectory_file', help='Binary trajectory file to plot')
+    parser.add_argument('--2d-only', action='store_true', dest='two_d_only', help='Only show 2D plot')
+    parser.add_argument('--3d-only', action='store_true', dest='three_d_only', help='Only show 3D plot')
+    parser.add_argument('--energy', action='store_true', help='Show energy plot')
+    parser.add_argument('--animate', action='store_true', help='Show animated trajectories')
+    parser.add_argument('--save-animation', type=str, help='Save animation as MP4 file')
+    parser.add_argument('--static', action='store_true', help='Show static plots (default if no --animate)')
+    parser.add_argument('--interval', type=int, help='Animation interval in milliseconds (default: auto-scaled to ~60s total)')
+    parser.add_argument('--target-duration', type=int, default=60, help='Target animation duration in seconds (default: 60)')
+    parser.add_argument('--center', type=str, help='Center animation on specified body (e.g., "Sun", "Earth")')
+    parser.add_argument('--list-bodies', action='store_true', help='List available bodies and exit')
+
+    args = parser.parse_args()
+    
+    print(f"Reading trajectory file: {args.trajectory_file}")
+    
+    try:
+        snapshots = read_trajectory_file(args.trajectory_file)
+        print(f"Loaded {len(snapshots)} snapshots")
+        
+        if not snapshots:
+            print("No data found in file!")
+            return
+        
+        trajectories, times = organize_trajectories(snapshots)
+        
+        # Handle list-bodies option
+        if args.list_bodies:
+            print(f"Available bodies in trajectory file:")
+            for body_name in sorted(trajectories.keys()):
+                print(f"  - {body_name}")
+            return
+        
+        print(f"Bodies found: {list(trajectories.keys())}")
+        print(f"Time range: {times[0]:.2e} - {times[-1]:.2e} seconds")
+        print(f"Number of time steps: {len(times)}")
+        
+        # Center trajectories on specified body if requested
+        original_trajectories = trajectories.copy()
+        if args.center:
+            try:
+                trajectories = center_trajectories_on_body(trajectories, args.center)
+                print(f"🎯 Centering animation on: {args.center}")
+            except ValueError as e:
+                print(f"❌ Error: {e}")
+                return
+        
+        # Check if we should animate or show static plots
+        show_animation = args.animate or args.save_animation
+        
+        if show_animation:
+            # Calculate animation parameters
+            interval_ms, anim_duration_sec, time_scale, is_manual = calculate_animation_params(
+                len(times), args.target_duration, args.interval
+            )
+            
+            print(f"\n🎬 Animation Settings:")
+            print(f"   Total frames: {len(times)}")
+            print(f"   Animation duration: {anim_duration_sec:.1f} seconds")
+            print(f"   Frame interval: {interval_ms}ms")
+            if is_manual:
+                print(f"   ⚙️  Using manual interval (--interval {args.interval})")
+                if time_scale != 1.0:
+                    print(f"   ⏱️  Time scale: {time_scale:.2f}x (animation {'faster' if time_scale > 1 else 'slower'} than target)")
+            else:
+                print(f"   🤖 Auto-scaled for {args.target_duration}s target duration")
+                print(f"   ⏱️  Time scale: 1.0x (optimized)")
+            
+            simulation_duration = times[-1] - times[0]
+            if simulation_duration > 0:
+                compression_ratio = simulation_duration / anim_duration_sec
+                if compression_ratio >= 86400:
+                    print(f"   📈 Compression: {compression_ratio/86400:.1f} days of simulation per second of animation")
+                elif compression_ratio >= 3600:
+                    print(f"   📈 Compression: {compression_ratio/3600:.1f} hours of simulation per second of animation")
+                else:
+                    print(f"   📈 Compression: {compression_ratio:.1f}x real-time")
+            print()
+            
+            print("Creating animated plots...")
+            animations = []
+            
+            # Create animated plots
+            if not args.three_d_only:
+                print("Creating 2D animation...")
+                fig_2d, anim_2d = plot_animated_2d(trajectories, times, interval_ms, args.center)
+                animations.append((fig_2d, anim_2d, '2d'))
+            
+            if not args.two_d_only:
+                print("Creating 3D animation...")
+                fig_3d, anim_3d = plot_animated_3d(trajectories, times, interval_ms, args.center)
+                animations.append((fig_3d, anim_3d, '3d'))
+            
+            # Save animations if requested
+            if args.save_animation:
+                for fig, anim, plot_type in animations:
+                    filename = f"{args.save_animation}_{plot_type}.mp4"
+                    print(f"Saving {plot_type.upper()} animation to {filename}...")
+                    try:
+                        anim.save(filename, writer='ffmpeg', fps=20)
+                        print(f"Animation saved to {filename}")
+                    except Exception as e:
+                        print(f"Error saving animation: {e}")
+                        print("Note: You may need to install ffmpeg for video export")
+            
+            plt.show()
+            
+        else:
+            print("Creating static plots...")
+            # Use original trajectories for static plots unless centering is requested
+            plot_trajectories = trajectories if args.center else original_trajectories
+            
+            # Plot static trajectories
+            if not args.three_d_only:
+                plot_trajectories_2d(plot_trajectories, times, args.center)
+            
+            if not args.two_d_only:
+                plot_trajectories_3d(plot_trajectories, times, args.center)
+            
+            if args.energy:
+                plot_energy_over_time(snapshots, times)
+            
+            plt.show()
+        
+    except FileNotFoundError:
+        print(f"Error: File '{args.trajectory_file}' not found!")
+    except Exception as e:
+        print(f"Error reading file: {e}")
+        import traceback
+        traceback.print_exc()
+
+if __name__ == "__main__":
+    main()

requirements.txt

@@ -0,0 +1,3 @@
+matplotlib>=3.5.0
+numpy>=1.20.0
+ffmpeg-python>=0.2.0

src/bin/orbiter.rs

@@ -14,7 +14,7 @@ struct CameraController {
     pub radius: f32,
     pub theta: f32, // azimuthal angle
     pub phi: f32,   // polar angle
-    pub last_mouse_pos: Option<Vec2>,
+    pub _last_mouse_pos: Option<Vec2>,
 }
 
 fn setup(
@@ -36,7 +36,7 @@ fn setup(
             radius,
             theta,
             phi,
-            last_mouse_pos: None,
+            _last_mouse_pos: None,
         },
     ));
     // Light
@@ -76,7 +76,7 @@ fn camera_orbit_controls(
     mouse_button_input: Res<ButtonInput<MouseButton>>,
     windows: Query<&Window>,
 ) {
-    let window = if let Some(window) = windows.iter().next() { window } else { return; };
+    let _window = if let Some(window) = windows.iter().next() { window } else { return; };
     let mut delta = Vec2::ZERO;
     for event in mouse_motion_events.read() {
         delta += event.delta;

src/bin/simulator.rs

@@ -1,6 +1,7 @@
 // Standard library
 use std::error::Error;
 use std::fs;
+use std::io::{self, Write};
 use std::path::Path;
 use std::time::{Duration,Instant};
 
@@ -9,10 +10,11 @@ use clap::Parser;
 use log::{info, debug};
 use indicatif::{ProgressBar, ProgressStyle};
 use humantime;
+use serde_json;
 
 // Internal modules from the library crate
 use orbital_simulator::simulation::Simulation;
-use orbital_simulator::types::norm_time;
+use orbital_simulator::types::{norm_time, real_time};
 use orbital_simulator as _; // for mod resolution
 
 #[derive(Parser, Debug)]
@@ -42,11 +44,27 @@ struct Args {
     /// Filename to save trajectory to
     #[arg(short, long)]
     output_file: String,
+
+    /// Force overwrite existing output file without confirmation
+    #[arg(short = 'w', long)]
+    force_overwrite: bool,
 }
 
 fn read_config<P: AsRef<Path>>(path: P) -> Result<orbital_simulator::config::ConfigFile, Box<dyn Error>> {
     let content = fs::read_to_string(&path)?;
-    let conf = toml::from_str(&content)?;
+    let path_str = path.as_ref().to_string_lossy();
+    
+    let conf: orbital_simulator::config::ConfigFile = if path_str.ends_with(".json") {
+        serde_json::from_str(&content)?
+    } else if path_str.ends_with(".toml") {
+        toml::from_str(&content)?
+    } else {
+        // Try JSON first, then TOML
+        serde_json::from_str(&content).or_else(|_| toml::from_str(&content))?
+    };
+
+    // Apply normalization settings if present
+    conf.apply_normalization();
 
     Ok(conf)
 }
@@ -73,6 +91,32 @@ fn format_duration_single_unit(dur: Duration) -> String {
     }
 }
 
+fn check_output_file(path: &str, force_overwrite: bool) -> Result<(), Box<dyn Error>> {
+    if Path::new(path).exists() {
+        if force_overwrite {
+            info!("Overwriting existing file: {}", path);
+            return Ok(());
+        }
+        
+        print!("⚠️  Output file '{}' already exists. Overwrite? [y/N]: ", path);
+        io::stdout().flush()?;
+        
+        let mut input = String::new();
+        io::stdin().read_line(&mut input)?;
+        
+        let response = input.trim().to_lowercase();
+        if response == "y" || response == "yes" {
+            info!("Overwriting existing file: {}", path);
+            Ok(())
+        } else {
+            println!("❌ Aborted. Use --force-overwrite (-w) to skip confirmation.");
+            std::process::exit(1);
+        }
+    } else {
+        Ok(())
+    }
+}
+
 fn main() {
     env_logger::init();
     let args = Args::parse();
@@ -80,13 +124,19 @@ fn main() {
     info!("Loading initial parameters from {}", args.config);
 
     let conf = read_config(args.config).unwrap();
+    
+    // Check if output file exists and handle overwrite
+    check_output_file(&args.output_file, args.force_overwrite).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);
     }
-
+    debug!("Simulation time and step size: {} seconds, {} seconds per step",
+        args.time.as_secs(), args.step_size);
     let n_steps = (args.time.as_secs() as f64 / args.step_size) as usize;
 
+    check_output_file(&args.output_file, args.force_overwrite).unwrap();
 
     let pb = ProgressBar::new(n_steps.try_into().unwrap());
     pb.set_style(
@@ -111,4 +161,4 @@ fn main() {
             pb.inc(args.steps_per_save as u64);
         })
     );
-} 
+}

src/config.rs

@@ -15,16 +15,24 @@ pub struct Body {
 
 #[derive(Debug, Deserialize, Serialize)]
 pub struct Normalization {
-    pub name: String,
+    pub m_0: f64,  // Mass normalization constant
-    pub mass: types::Mass,
+    pub t_0: f64,  // Time normalization constant  
-    pub position: types::Position,
+    pub r_0: f64,  // Distance normalization constant
-    pub velocity: types::Velocity,
-
-    #[serde(default)]
-    pub acceleration: types::Acceleration,
 }
 
 #[derive(Debug, Deserialize)]
 pub struct ConfigFile {
     pub bodies: Vec<Body>,
+    
+    #[serde(default)]
+    pub normalization: Option<Normalization>,
+}
+
+impl ConfigFile {
+    /// Apply normalization settings if present in the config
+    pub fn apply_normalization(&self) {
+        if let Some(ref norm) = self.normalization {
+            types::set_normalization(norm);
+        }
+    }
 }

src/simulation.rs

@@ -34,7 +34,7 @@ impl Simulation {
             bodies[i].position = norm_pos(bodies[i].position);
             bodies[i].velocity = norm_vel(bodies[i].velocity);
         }
-        step_size = norm_time(step_size);
+        step_size = step_size;
         Self {bodies, step_size, steps_per_save, save_file}
     }
 
@@ -42,12 +42,15 @@ impl Simulation {
     where F: FnMut() {
         let file = OpenOptions::new()
             .create(true)
-            .append(true)
+            .write(true)
+            .truncate(true)
             .open(&self.save_file);
         let mut writer = BufWriter::new(file.unwrap());
 
         for i in 0..steps {
-            debug!("Step: {}", i);
+            trace!("norm step size is {}", self.step_size);
+            trace!("real step size is {}", real_time(self.step_size));
+            trace!("Step: {}, Time: {}", i, real_time(i as f64 * self.step_size));
             self.reset_accelerations();
             self.calculate_accelerations();
             self.step_bodies();
@@ -55,7 +58,7 @@ impl Simulation {
                 //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),
+                    time: real_time(i as f64 * self.step_size),
                     bodies: &real_bodies,
                 };
 
@@ -75,8 +78,9 @@ impl Simulation {
             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;
+                // In normalized units, G = 1
-                self.bodies[j].acceleration -= r_hat_over_r3[i][j]*mass_i;
+                self.bodies[i].acceleration += r_hat_over_r3[i][j] * mass_j;
+                self.bodies[j].acceleration -= r_hat_over_r3[i][j] * mass_i;
             }
         }
     }

src/types.rs

@@ -1,8 +1,10 @@
 #![allow(unused)]
 use glam::DVec3;
 use once_cell::sync::Lazy;
+use std::sync::RwLock;
+use log::{debug};
 
-use crate::config::Body;
+use crate::config::{Body, Normalization};
 
 pub type Position = DVec3;
 pub type Velocity = DVec3;
@@ -24,62 +26,119 @@ 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(|| {
+// Default normalization constants
-    (R_0.powf(3.0) / (G * M_0)).sqrt()
+const DEFAULT_R_0: f64 = EARTH_RADIUS;
+const DEFAULT_M_0: f64 = EARTH_MASS;
+
+// Global normalization context
+static NORMALIZATION: Lazy<RwLock<NormalizationContext>> = Lazy::new(|| {
+    RwLock::new(NormalizationContext::default())
 });
 
+#[derive(Debug, Clone)]
+pub struct NormalizationContext {
+    pub r_0: f64,
+    pub m_0: f64,
+    pub t_0: f64,
+}
+
+impl Default for NormalizationContext {
+    fn default() -> Self {
+        let r_0 = DEFAULT_R_0;
+        let m_0 = DEFAULT_M_0;
+        let t_0 = (r_0.powf(3.0) / (G * m_0)).sqrt();
+        debug!("Using default normalization: r_0 = {}, m_0 = {}, t_0 = {}", r_0, m_0, t_0);
+        Self { r_0, m_0, t_0 }
+    }
+}
+
+impl From<&Normalization> for NormalizationContext {
+    fn from(norm: &Normalization) -> Self {
+        Self { 
+            r_0: norm.r_0, 
+            m_0: norm.m_0, 
+            t_0: norm.t_0 
+        }
+    }
+}
+
+/// Set the global normalization context from a config normalization
+pub fn set_normalization(norm: &Normalization) {
+    let mut context = NORMALIZATION.write().unwrap();
+    *context = NormalizationContext::from(norm);
+}
+
+/// Reset to default normalization
+pub fn reset_normalization() {
+    let mut context = NORMALIZATION.write().unwrap();
+    *context = NormalizationContext::default();
+}
+
+/// Get current normalization context
+pub fn get_normalization() -> NormalizationContext {
+    NORMALIZATION.read().unwrap().clone()
+}
+
 
 #[inline]
 pub fn norm_pos(pos: Position) -> Position {
-    pos / R_0
+    let norm = get_normalization();
+    pos / norm.r_0
 }
 
 #[inline]
 pub fn real_pos(pos: Position) -> Position {
-    pos * R_0
+    let norm = get_normalization();
+    pos * norm.r_0
 }
 
 #[inline]
 pub fn norm_mass(mass: Mass) -> Mass {
-    mass / M_0
+    let norm = get_normalization();
+    mass / norm.m_0
 }
 
 #[inline]
 pub fn real_mass(mass: Mass) -> Mass {
-    mass * M_0
+    let norm = get_normalization();
+    mass * norm.m_0
 }
 
 #[inline]
 pub fn norm_time(time: Time) -> Time {
-    time / *T_0
+    let norm = get_normalization();
+    time / norm.t_0
 }
 
 #[inline]
 pub fn real_time(time: Time) -> Time {
-    time * *T_0
+    let norm = get_normalization();
+    time * norm.t_0
 }
 
 #[inline]
 pub fn norm_vel(vel: Velocity) -> Velocity {
-    vel / (R_0 / *T_0)
+    let norm = get_normalization();
+    vel / (norm.r_0 / norm.t_0)
 }
 
 #[inline]
 pub fn real_vel(vel: Velocity) -> Velocity {
-    vel * (R_0 / *T_0)
+    let norm = get_normalization();
+    vel * (norm.r_0 / norm.t_0)
 }
 
 #[inline]
 pub fn norm_acc(acc: Acceleration) -> Acceleration {
-    acc / (R_0 / (*T_0 * *T_0))
+    let norm = get_normalization();
+    acc / (norm.r_0 / (norm.t_0 * norm.t_0))
 }
 
 #[inline]
 pub fn real_acc(acc: Acceleration) -> Acceleration {
-    acc * (R_0 / (*T_0 * *T_0))
+    let norm = get_normalization();
+    acc * (norm.r_0 / (norm.t_0 * norm.t_0))
 }
 
 #[inline]