diff --git a/orbiter/orbits/body.py b/orbiter/orbits/body.py index 57028a6..67189d6 100644 --- a/orbiter/orbits/body.py +++ b/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}") + diff --git a/orbiter/orbits/simulator.py b/orbiter/orbits/simulator.py index 76337dc..ee1fb6d 100644 --- a/orbiter/orbits/simulator.py +++ b/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() - self.move_bodies() - if (steps % self.steps_per_save == 0): - self._checkpoint() + 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() 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 = G*self.bodies[i].m*self.bodies[j].m/(dists[i][j]**2) - print(vec_norm, f) - force_sum += vec_norm*f + f = self.bodies[i].m*self.bodies[j].m/(dists[i][j]**2) + force_sum += vec_norm*f self.bodies[i].A = force_sum/self.bodies[i].m def move_bodies(self): diff --git a/orbiter/units.py b/orbiter/units.py index 9ab7ee7..3a37d5a 100644 --- a/orbiter/units.py +++ b/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) + + diff --git a/test.py b/test.py new file mode 100644 index 0000000..922f52f --- /dev/null +++ b/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)) \ No newline at end of file