tried converting to Decimal for high precision

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

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
 

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)
 
 

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