Having magnitude issues - of course!

orbiter/orbits/body.py

@@ -7,7 +7,7 @@ class Body:
     """
     Base class for any orbital body
     """
-    def __init__(self, X: Position, V: Velocity, m: Mass):
+    def __init__(self, X: Position, V: Velocity, m: Mass, name: str = ""):
         """
         x (position) and v (velocity) 
         """
@@ -15,6 +15,7 @@ class Body:
         self.V = V
         self.A = Acceleration([0,0,0])
         self.m = m
+        self.name = name
 
     def save(self):
         return (self.X, self.V, self.m)
@@ -27,4 +28,7 @@ class Body:
         self.X = step_size*self.V
         self.V = step_size*self.A
     
+    def __str__(self):
+        return str(f"{self.name}: X = {self.X}, V = {self.V}")
+    
 

orbiter/orbits/simulator.py

@@ -4,7 +4,7 @@ from scipy.spatial import distance_matrix
 
 from .body import Body
 
-G = 1
+
 
 class Simulator:
     """
@@ -104,10 +104,13 @@ class Simulator:
                 steps=stepsz_n_np)
 
     def run(self, steps):
-        self.calculate_forces()
+        for i in range(steps):
-        self.move_bodies()
+            print(f"On the {i}th step - vel {self.bodies[1].V}")
-        if (steps % self.steps_per_save == 0):
+            self.calculate_forces()
-            self._checkpoint()
+            self.move_bodies()
+            self.current_step += 1
+            if (self.current_step % self.steps_per_save == 0):
+                self._checkpoint()
 
     def calculate_forces(self):
         positions = [
@@ -115,16 +118,16 @@ class Simulator:
         ]
         dists = distance_matrix(positions, positions)
         for i in range(len(self.bodies)):
-            force_sum = np.array([0.0,0.0,0.0])
+            print(str(self.bodies[i]))
+            force_sum = np.array([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
-                print(vec)
                 norm = np.linalg.norm(vec)
                 vec_norm = vec/norm
-                if (i != j):
+                f = self.bodies[i].m*self.bodies[j].m/(dists[i][j]**2)
-                    f = G*self.bodies[i].m*self.bodies[j].m/(dists[i][j]**2)
+                force_sum += vec_norm*f
-                    print(vec_norm, f)
-                    force_sum += vec_norm*f
             self.bodies[i].A = force_sum/self.bodies[i].m
 
     def move_bodies(self):

orbiter/units.py

@@ -1,9 +1,53 @@
 import numpy as np
 
+#Declare these as units to make code clearer
 Position = np.array
-
 Velocity = np.array
-
 Acceleration = np.array
-
 Mass = int
+
+####################
+# USEFUL CONSTANTS
+####################
+EARTH_MASS = 5.972e24 #kg
+EARTH_RADIUS = 6.378e6 #meters
+EARTH_ORBITAL_VELOCITY = 29.78e3
+AU = 149_597_870_700 #meters
+
+MOON_MASS = 7.34767309e22
+MOON_ORBITAL_VELOCITY = 1.022e3 #m/s relative to earth
+
+SUN_MASS = 1.989e30 #kg
+SUN_RADIUS = 6.957e8 #meters
+
+#NORMALIZING CONSTANTS
+G = 6.67430e-11
+r_0 = EARTH_RADIUS #1.496e11
+m_0 = 5.972e24
+t_0 = np.sqrt((r_0**3) / (G*m_0))
+
+def norm_pos(pos):
+    return pos / r_0
+
+def real_pos(pos):
+    return pos * r_0
+
+def norm_mass(mass):
+    return mass / m_0
+
+def real_mass(mass):
+    return mass * m_0
+
+def norm_time(time):
+    return time / t_0
+
+def real_time(time):
+    return time * t_0
+
+def norm_vel(vel):
+    return vel / (r_0/t_0)
+
+def real_vel(vel):
+    return vel * (r_0/t_0)
+
+

test.py

@@ -0,0 +1,31 @@
+from orbiter.orbits.body import Body
+from orbiter.orbits.simulator import Simulator
+from orbiter.units import *
+
+from pathlib import Path
+
+#set up the earth
+earth = Body(
+    Position([0,0,0]), 
+    Velocity([0,0,0]), 
+    Mass(norm_mass(EARTH_MASS)),
+    "Earth"
+)
+
+#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
+    Mass(norm_mass(80)), #avg person
+    "Person"
+)
+
+time_to_run = norm_time(5)
+STEP_SIZE = 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))
+
+print(real_vel(person.V))