Commit 498165e1 authored by Olaf's avatar Olaf

Add option to not save data at initial time

parent d4fa1921
This diff is collapsed.
......@@ -84,7 +84,7 @@ def calc(bodies, settings):
pass
else:
pass
positions, velocities, times, E_kin = dynamics(bodies, settings.dimensions, settings.timestep, settings.time, settings.barneshut, settings.theta, settings.approx, settings.timesave)
positions, velocities, times, E_kin = dynamics(bodies, settings.dimensions, settings.timestep, settings.time, settings.barneshut, settings.theta, settings.approx, settings.timesave, settings.startsave)
for i, body in enumerate(bodies):
if isinstance(body, Body):
body.add_trajectory((positions.transpose(1,0,2)[i]).tolist())
......
......@@ -225,7 +225,7 @@ def force_approx(mass, positions, D, n):
return Force
def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave):
"""
Calculates the dynamics (trajectories) of all bodies
......@@ -247,6 +247,8 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
states if only interaction on asteroids and not of the asteroids should be taken
j: int
every j_th timestep the data is saved
startsave: int
start saving at this timestep
Returns:
--------
......@@ -262,6 +264,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
T = round(t_max/h)
positions, velo, mass, n = extract_data(bodies, D)
if barnes_hut and theta != 0 and not approximation:
force0 = force_barneshut(positions, mass, theta)
elif approximation and not barnes_hut:
......@@ -269,11 +272,16 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
else:
force0 = force_calc(mass, positions, D)
all_pos = np.zeros((math.ceil((T+1)/j), positions.shape[0], positions.shape[1]))
times = [0]
all_pos = np.zeros((math.ceil((T-startsave+1)/j), positions.shape[0], positions.shape[1]))
E_kin = []
all_pos[0] = positions
index = 0
times = []
index = -1
if startsave == 0:
index += 1
times.append(0)
E_kin.append(0.5*mass*(np.linalg.norm(velo, axis=1)**2))
all_pos[0] = positions
if barnes_hut and theta != 0 and not approximation:
for i in range(1,T+1):
......@@ -281,7 +289,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
force1 = force_barneshut(positions, mass, theta)
velo = velocity(velo, force0, force1, h)
force0 = force1
if i%j == 0:
if i%j == 0 and i>=startsave:
index += 1
all_pos[index] = positions
E_kin.append(0.5*mass*(np.linalg.norm(velo, axis=1)**2))
......@@ -292,7 +300,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
force1 = force_approx(mass, positions, D, n)
velo = velocity(velo, force0, force1, h)
force0 = force1
if i%j == 0:
if i%j == 0 and i>=startsave:
index += 1
all_pos[index] = positions
E_kin.append(0.5*mass*(np.linalg.norm(velo, axis=1)**2))
......@@ -303,7 +311,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j):
force1 = force_calc(mass, positions, D)
velo = velocity(velo, force0, force1, h)
force0 = force1
if i%j == 0:
if i%j == 0 and i>=startsave:
index += 1
all_pos[index] = positions
E_kin.append(0.5*mass*(np.linalg.norm(velo, axis=1)**2))
......
......@@ -522,6 +522,11 @@ class Settings:
self.timesave = sets['save_timestep']
except:
self.timesave = 1
try:
self.startsave = sets['start_save']
except:
self.startsave = 0
......
......@@ -150,7 +150,7 @@ Neptune = Body(Neptune_parameters); solar_system.append(Neptune)
# Asteroid belt, data from: http://astronomy.swin.edu.au/~gmackie/SAO/HET602_Amaya/m14a01.ppt
Asteroid_belt_parameters = {
'N' : 10, # number of asteroids
'N' : 100, # number of asteroids
'mass' : 0.012, # total mass of all asteroids combined
'R_min' : 2.1, # lower bound orbital radius from the Sun
'R_max' : 4.1, # upper bound orbital radius from the Sun
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment