it WORKS!

orbiter/orbits/body.py

@@ -1,7 +1,8 @@
 from enum import Enum
 import numpy as np
 
-from ..units import Position, Velocity, Mass, Acceleration
+from ..units import *
+from .calc import format_sig_figs
 
 class Body:
     """
@@ -25,10 +26,32 @@ class Body:
         return cls(tup[0], tup[1], tup[2])
 
     def step(self, step_size: float):
-        self.X = step_size*self.V
+        self.X += step_size*self.V
-        self.V = step_size*self.A
+        self.V += step_size*self.A
+        self.A = HighPrecisionVector([0,0,0])
     
     def __str__(self):
-        return str(f"{self.name}: X = {self.X}, V = {self.V}")
+        pos = 'm '.join([format_sig_figs(real_pos(x), 3) for x in self.X])
+        vel = 'm/s '.join([format_sig_figs(real_vel(v), 3) for v in self.V])
+        return str(f"{self.name}: X = {pos}, V = {vel}")
+    
+    def E(self):
+        return self.ke() + self.pe()
+
+    def pe(self):
+        return -self.m/self.dist_from_o()
+
+    def dist_from_o(self):
+        return sum([x**2 for x in self.X]).sqrt()
+
+
+    def ke(self):
+        return Decimal(0.5)*self.m*(self._speed()**2)
+
+    def _speed(self):
+        return sum([v**2 for v in self.V]).sqrt()
+
+    def speed(self):
+        return str(f"{format_sig_figs(real_vel(self._speed),5)}")
     
 

orbiter/orbits/calc.py

@@ -1,5 +1,7 @@
 import numpy as np
+import sys
 from decimal import Decimal
+from time import time
 
 from ..units import HighPrecisionMatrix
 
@@ -13,3 +15,63 @@ def calculate_distances(positions):
             dists[i][j] = Decimal(d)
             dists[j][i] = Decimal(d)
     return dists
+
+def print_progress_bar(iteration, total, start_time, length=50):
+    """Prints a progress bar to the console."""
+    percent = (iteration / total) * 100
+    filled_length = int(length * iteration // total)
+    bar = '#' * filled_length + '-' * (length - filled_length)
+    sys.stdout.write(f'\r[{bar}] {percent:.2f}% {int(iteration/(time()-start_time))} steps/s')
+    sys.stdout.flush()
+
+
+def format_sig_figs(value, sig_figs):
+    """Format a number to a specified number of significant figures for printing."""
+    if value == 0:
+        return "0"
+    return f"{value:.{sig_figs-1}e}"
+
+
+
+def plot_points_terminal(vectors, stdscr, scale=500000, grid_width=30, grid_height=30):
+    """Plots multiple points in the terminal, scaled and centered at (0,0)."""
+    stdscr.clear()
+
+    if not vectors:
+        stdscr.addstr(0, 0, "No vectors provided.")
+        stdscr.refresh()
+        return
+
+    # Apply scaling
+    scaled_vectors = {(round(vec[0] / scale), round(vec[1] / scale)) for vec in vectors}
+
+    # Find min and max coordinates
+    min_x = min(vec[0] for vec in scaled_vectors)
+    max_x = max(vec[0] for vec in scaled_vectors)
+    min_y = min(vec[1] for vec in scaled_vectors)
+    max_y = max(vec[1] for vec in scaled_vectors)
+
+    # Center offsets to keep (0,0) in middle
+    center_x = (grid_width // 2) - min_x
+    center_y = (grid_height // 2) - min_y
+
+    # Adjust coordinates for plotting
+    adjusted_vectors = {(vec[0] + center_x, vec[1] + center_y) for vec in scaled_vectors}
+
+    # Ensure grid boundaries
+    max_terminal_y, max_terminal_x = stdscr.getmaxyx()
+    max_x = min(grid_width, max_terminal_x - 5)
+    max_y = min(grid_height, max_terminal_y - 5)
+
+    # Draw grid with points
+    for i in range(grid_height, -1, -1):
+        row = f"{i - center_y:2} | "
+        for j in range(grid_width + 1):
+            row += "● " if (j, i) in adjusted_vectors else ". "
+        stdscr.addstr(max_y - i, 0, row[:max_terminal_x - 1])  # Ensure no overflow
+
+    # Print X-axis labels
+    x_labels = "   " + " ".join(f"{j - center_x:2}" for j in range(max_x + 1))
+    stdscr.addstr(max_y + 1, 0, x_labels[:max_terminal_x - 1])  # Avoid out-of-bounds error
+
+    stdscr.refresh()

orbiter/orbits/simulator.py

@@ -1,10 +1,11 @@
 from pathlib import Path
 import numpy as np
 from decimal import InvalidOperation, DivisionByZero
+from time import time
 
 from .body import Body
-from .calc import calculate_distances
+from .calc import *
-from ..units import HighPrecisionVector
+from ..units import *
 
 
 
@@ -106,13 +107,31 @@ class Simulator:
                 steps=stepsz_n_np)
 
     def run(self, steps):
+        time_start = time()
+        import matplotlib.pyplot as plt
+        plt.ion()
+        fig, ax = plt.subplots()
+        x_data, y_data =[],[]
+        line, = ax.plot(x_data, y_data, 'bo-')
+        ax.set_xlim(-8e6,8e6)
+        ax.set_ylim(-8e6,8e6)
         for i in range(steps):
-            #print(f"On the {i}th step - vel {self.bodies[1].V}")
             self.calculate_forces()
             self.move_bodies()
             self.current_step += 1
             if (self.current_step % self.steps_per_save == 0):
-                self._checkpoint()
+                for b in self.bodies:
+                    x_data.append(int(real_pos(b.X[0])))
+                    y_data.append(int(real_pos(b.X[1]))) 
+                line.set_xdata(x_data)
+                line.set_ydata(y_data)
+                plt.draw()
+                plt.pause(0.2)
+                x_data, y_data =[],[]
+
+
+                #print_progress_bar(i, steps, time_start)
+                #self._checkpoint()
 
     def calculate_forces(self):
         positions = [
@@ -120,21 +139,14 @@ class Simulator:
         ]
         dists = calculate_distances(positions)
         for i in range(len(self.bodies)):
-            print(str(self.bodies[i]))
+            for j in range(i, len(self.bodies)):
-            force_sum = HighPrecisionVector([0.0, 0.0, 0.0])
-            for j in range(len(self.bodies)):
                 if i == j:
                     continue
                 vec = self.bodies[i].X - self.bodies[j].X
-                norm = np.linalg.norm(vec)
+                f = vec/(dists[i][j]**3)
-                try:
+                self.bodies[i].A += f*self.bodies[j].m
-                    vec_norm = vec/norm
+                self.bodies[j].A += f*self.bodies[i].m
-                    f = self.bodies[i].m*self.bodies[j].m/(dists[i][j]**2)
+
-                except (InvalidOperation, DivisionByZero):
-                    vec_norm = HighPrecisionVector([0,0,0])
-                    f = 0
-                force_sum += vec_norm*f
-            self.bodies[i].A = force_sum/self.bodies[i].m
 
     def move_bodies(self):
         for body in self.bodies:

orbiter/units.py

@@ -10,7 +10,6 @@ class HighPrecisionVector(np.ndarray):
 
 class HighPrecisionMatrix(np.ndarray):
     def __new__(cls, dim1, dim2, *args, **kwargs):
-        print(f"dim1{dim1}, dim2{dim2}")
         decimal_array = [Decimal(0) for _ in range(dim1)]
         decimal_matrix = [decimal_array for _ in range(dim2)]
         obj = np.asarray(decimal_matrix).view(cls)
@@ -39,6 +38,9 @@ MOON_ORBITAL_VELOCITY = Decimal(1022) #m/s relative to earth
 SUN_MASS = Decimal(1989 * 10**27) #kg
 SUN_RADIUS = Decimal(6957 * 10**5) #meters
 
+pi_approx = Decimal("3.14159265358979323846264338327950288419716939937510")
+
+
 #NORMALIZING CONSTANTS
 G = Decimal(6.67430e-11)
 r_0 = Decimal(EARTH_RADIUS) #1.496e11
@@ -49,7 +51,7 @@ def norm_pos(pos):
     return Decimal(pos) / Decimal(r_0)
 
 def real_pos(pos):
-    return Decimal(pos) * Decimal(r_0)
+    return pos * Decimal(r_0)
 
 def norm_mass(mass):
     return Decimal(mass) / Decimal(m_0)
@@ -64,9 +66,9 @@ def real_time(time):
     return Decimal(time) * Decimal(t_0)
 
 def norm_vel(vel):
-    return Decimal(vel) / Decimal(r_0/t_0)
+    return vel / Decimal(r_0/t_0)
 
 def real_vel(vel):
-    return Decimal(vel) * Decimal(r_0/t_0)
+    return vel * Decimal(r_0/t_0)
 
 

test.py

@@ -5,7 +5,7 @@ from orbiter.units import *
 from pathlib import Path
 from decimal import Decimal, getcontext
 
-getcontext().prec = 100
+getcontext().prec = 50
 
 #set up the earth
 earth = Body(
@@ -15,19 +15,29 @@ earth = Body(
     "Earth"
 )
 
+r = EARTH_RADIUS+100_000
 #Lets try a body just outside earth accelerating in. Should be 9.8m/s2
 person = Body(
-    Position([norm_pos(EARTH_RADIUS+100_000),0,0]), #10_000m in the sky, airliner height!
+    Position([norm_pos(r),0,0]), #10_000m in the sky, airliner height!
-    Velocity([0,0,0]), #start from standstill
+    Velocity([0,(Decimal(0.5)/norm_pos(r)).sqrt(),0]), #orbital velocity
     Mass(norm_mass(80)), #avg person
     "Person"
 )
 
-time_to_run = norm_time(5)
+T = 2*pi_approx*norm_pos(r)/person.V[1]
-STEP_SIZE = Decimal(1e-10)
-n_steps = time_to_run/STEP_SIZE
 
-s = Simulator([earth,person], STEP_SIZE, 1, Path("hello_world"))
+time_to_run = 15 #norm_time(2000)
-s.run(10)
+STEP_SIZE = Decimal(6e-4)
+n_steps = int(time_to_run/STEP_SIZE)
 
-print(real_vel(person.V))
+def main():
+    print("Before: ")
+    print(str(person))
+    print(str(earth))
+    s = Simulator([earth,person], STEP_SIZE, 100, Path("hello_world"))
+    s.run(n_steps)
+    print("\nAfter:")
+    print(str(person))
+    print(str(earth))
+
+main()