From cf81cd17e61a9e72a5689c1e97d1e8933840dc40 Mon Sep 17 00:00:00 2001 From: Thomas Faour Date: Sat, 31 May 2025 17:38:26 -0400 Subject: [PATCH] tried converting to Decimal for high precision --- orbiter/orbits/calc.py | 15 +++++++-- orbiter/orbits/simulator.py | 18 ++++++---- orbiter/units.py | 65 ++++++++++++++++++++++++------------- test.py | 8 +++-- 4 files changed, 72 insertions(+), 34 deletions(-) diff --git a/orbiter/orbits/calc.py b/orbiter/orbits/calc.py index 4558b02..51ddfb7 100644 --- a/orbiter/orbits/calc.py +++ b/orbiter/orbits/calc.py @@ -1,4 +1,15 @@ -from scipy.spatial import distance_matrix +import numpy as np +from decimal import Decimal + +from ..units import HighPrecisionMatrix def calculate_distances(positions): - return distance_matrix(positions, positions) \ No newline at end of file + N = len(positions) + dists = HighPrecisionMatrix(N,N) + for i in range(N): + dists[i][i] = Decimal(0) + for j in range(i+1, N): + d = sum([x**2 for x in (positions[i]-positions[j])]).sqrt() + dists[i][j] = Decimal(d) + dists[j][i] = Decimal(d) + return dists diff --git a/orbiter/orbits/simulator.py b/orbiter/orbits/simulator.py index ee1fb6d..07918f6 100644 --- a/orbiter/orbits/simulator.py +++ b/orbiter/orbits/simulator.py @@ -1,8 +1,10 @@ from pathlib import Path import numpy as np -from scipy.spatial import distance_matrix +from decimal import InvalidOperation, DivisionByZero from .body import Body +from .calc import calculate_distances +from ..units import HighPrecisionVector @@ -105,7 +107,7 @@ class Simulator: def run(self, steps): for i in range(steps): - print(f"On the {i}th step - vel {self.bodies[1].V}") + #print(f"On the {i}th step - vel {self.bodies[1].V}") self.calculate_forces() self.move_bodies() self.current_step += 1 @@ -116,17 +118,21 @@ class Simulator: positions = [ body.X for body in self.bodies ] - dists = distance_matrix(positions, positions) + dists = calculate_distances(positions) for i in range(len(self.bodies)): print(str(self.bodies[i])) - force_sum = np.array([0.0, 0.0, 0.0]) + 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) - vec_norm = vec/norm - f = self.bodies[i].m*self.bodies[j].m/(dists[i][j]**2) + 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 diff --git a/orbiter/units.py b/orbiter/units.py index 3a37d5a..e9eab20 100644 --- a/orbiter/units.py +++ b/orbiter/units.py @@ -1,53 +1,72 @@ import numpy as np +from decimal import Decimal #Declare these as units to make code clearer -Position = np.array -Velocity = np.array -Acceleration = np.array -Mass = int +class HighPrecisionVector(np.ndarray): + def __new__(cls, input_array, *args, **kwargs): + decimal_array = [Decimal(i) for i in input_array] + obj = np.asarray(decimal_array).view(cls) + return obj + +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) + return obj + +class Position(HighPrecisionVector): + pass +class Velocity(HighPrecisionVector): + pass +class Acceleration(HighPrecisionVector): + pass + +Mass = Decimal #################### # USEFUL CONSTANTS #################### -EARTH_MASS = 5.972e24 #kg -EARTH_RADIUS = 6.378e6 #meters -EARTH_ORBITAL_VELOCITY = 29.78e3 -AU = 149_597_870_700 #meters +EARTH_MASS = Decimal(5972 * 10**21) #kg +EARTH_RADIUS = Decimal(6378 * 10**3) #meters +EARTH_ORBITAL_VELOCITY = Decimal(29780) # m/s +AU = Decimal(149_597_870_700) #meters -MOON_MASS = 7.34767309e22 -MOON_ORBITAL_VELOCITY = 1.022e3 #m/s relative to earth +MOON_MASS = Decimal(734767309 * 10**14) +MOON_ORBITAL_VELOCITY = Decimal(1022) #m/s relative to earth -SUN_MASS = 1.989e30 #kg -SUN_RADIUS = 6.957e8 #meters +SUN_MASS = Decimal(1989 * 10**27) #kg +SUN_RADIUS = Decimal(6957 * 10**5) #meters #NORMALIZING CONSTANTS -G = 6.67430e-11 -r_0 = EARTH_RADIUS #1.496e11 -m_0 = 5.972e24 +G = Decimal(6.67430e-11) +r_0 = Decimal(EARTH_RADIUS) #1.496e11 +m_0 = Decimal(5.972e24) t_0 = np.sqrt((r_0**3) / (G*m_0)) def norm_pos(pos): - return pos / r_0 + return Decimal(pos) / Decimal(r_0) def real_pos(pos): - return pos * r_0 + return Decimal(pos) * Decimal(r_0) def norm_mass(mass): - return mass / m_0 + return Decimal(mass) / Decimal(m_0) def real_mass(mass): - return mass * m_0 + return Decimal(mass) * Decimal(m_0) def norm_time(time): - return time / t_0 + return Decimal(time) / Decimal(t_0) def real_time(time): - return time * t_0 + return Decimal(time) * Decimal(t_0) def norm_vel(vel): - return vel / (r_0/t_0) + return Decimal(vel) / Decimal(r_0/t_0) def real_vel(vel): - return vel * (r_0/t_0) + return Decimal(vel) * Decimal(r_0/t_0) diff --git a/test.py b/test.py index 922f52f..ca4df28 100644 --- a/test.py +++ b/test.py @@ -3,6 +3,9 @@ from orbiter.orbits.simulator import Simulator from orbiter.units import * from pathlib import Path +from decimal import Decimal, getcontext + +getcontext().prec = 100 #set up the earth earth = Body( @@ -21,11 +24,10 @@ person = Body( ) time_to_run = norm_time(5) -STEP_SIZE = 1e-10 +STEP_SIZE = Decimal(1e-10) n_steps = time_to_run/STEP_SIZE s = Simulator([earth,person], STEP_SIZE, 1, Path("hello_world")) -print(n_steps) -s.run(int(time_to_run/STEP_SIZE)) +s.run(10) print(real_vel(person.V)) \ No newline at end of file