From 5e3d043af4e2419323adb9c68a06205e14c822ff Mon Sep 17 00:00:00 2001 From: Thomas Faour Date: Sun, 1 Jun 2025 00:42:31 -0400 Subject: [PATCH] it WORKS! --- orbiter/orbits/body.py | 31 ++++++++++++++++--- orbiter/orbits/calc.py | 62 +++++++++++++++++++++++++++++++++++++ orbiter/orbits/simulator.py | 44 ++++++++++++++++---------- orbiter/units.py | 10 +++--- test.py | 28 +++++++++++------ 5 files changed, 142 insertions(+), 33 deletions(-) diff --git a/orbiter/orbits/body.py b/orbiter/orbits/body.py index 67189d6..676245f 100644 --- a/orbiter/orbits/body.py +++ b/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.V = step_size*self.A + self.X += step_size*self.V + 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)}") diff --git a/orbiter/orbits/calc.py b/orbiter/orbits/calc.py index 51ddfb7..71c0e8a 100644 --- a/orbiter/orbits/calc.py +++ b/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() diff --git a/orbiter/orbits/simulator.py b/orbiter/orbits/simulator.py index 07918f6..4ae4362 100644 --- a/orbiter/orbits/simulator.py +++ b/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 ..units import HighPrecisionVector +from .calc import * +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])) - force_sum = HighPrecisionVector([0.0, 0.0, 0.0]) - for j in range(len(self.bodies)): + for j in range(i, len(self.bodies)): if i == j: continue vec = self.bodies[i].X - self.bodies[j].X - norm = np.linalg.norm(vec) - try: - vec_norm = vec/norm - 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 + f = vec/(dists[i][j]**3) + self.bodies[i].A += f*self.bodies[j].m + self.bodies[j].A += f*self.bodies[i].m + def move_bodies(self): for body in self.bodies: diff --git a/orbiter/units.py b/orbiter/units.py index e9eab20..be121ae 100644 --- a/orbiter/units.py +++ b/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) diff --git a/test.py b/test.py index ca4df28..cc8ed07 100644 --- a/test.py +++ b/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! - Velocity([0,0,0]), #start from standstill + Position([norm_pos(r),0,0]), #10_000m in the sky, airliner height! + Velocity([0,(Decimal(0.5)/norm_pos(r)).sqrt(),0]), #orbital velocity Mass(norm_mass(80)), #avg person "Person" ) -time_to_run = norm_time(5) -STEP_SIZE = Decimal(1e-10) -n_steps = time_to_run/STEP_SIZE +T = 2*pi_approx*norm_pos(r)/person.V[1] -s = Simulator([earth,person], STEP_SIZE, 1, Path("hello_world")) -s.run(10) +time_to_run = 15 #norm_time(2000) +STEP_SIZE = Decimal(6e-4) +n_steps = int(time_to_run/STEP_SIZE) -print(real_vel(person.V)) \ No newline at end of file +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() \ No newline at end of file