Commit 283d1b9f authored by Olaf's avatar Olaf

Now able to do calculations based on previous data

parent c4444803
This diff is collapsed.
......@@ -86,7 +86,8 @@ def calc(bodies, settings):
pass
positions, velocities, times, E_kin, E_pot = dynamics(bodies, settings.dimensions, settings.timestep,
settings.time, settings.barneshut, settings.theta,
settings.approx, settings.timesave, settings.startsave)
settings.approx, settings.timesave, settings.startsave,
settings.continuecalc, settings.filename)
for i, body in enumerate(bodies):
if isinstance(body, Body):
body.add_trajectory((positions.transpose(1,0,2)[i]).tolist())
......
......@@ -3,7 +3,7 @@ import warnings
from predef import *
from barnes_hut import *
import math
import matplotlib.pyplot as plt
import pickle
def extract_data(bodies, D):
......@@ -85,6 +85,63 @@ def extract_data(bodies, D):
# return
return positions, velocities, masses, n
def extract_file_data(file):
"""
Extracts data from a file to continue calculations.
Returns:
pos: all positions
velo: all velocities
mass: all masses
n: number of non-asteroid bodies
t: final time of last calculation
"""
# Check whether file extension is given
if file[-4:] != ".pkl":
file += ".pkl"
with open(file, 'rb') as f:
data = pickle.load(f)
f.close()
bodies = data['bodies']
settings_saved = data['settings']
barnes_hut = settings_saved.barneshut
theta = settings_saved.theta
approximation = settings_saved.approx
t = bodies[0].times[-1]
velo = []
pos = []
mass = []
n = 0
for i in range(len(bodies)):
if isinstance(bodies[i], Asteroids):
if bodies[i].N >1:
traj = np.asarray(bodies[i].trajectories)
traj = np.transpose(traj, (1,0,2))
pos.extend((traj[-1]).tolist())
velo.extend(bodies[i].velocities)
mass.extend(bodies[i].masses)
elif bodies[i].N == 1:
pos.append(bodies[i].trajectories[-1])
velo.append(bodies[i].velocities)
mass.append(bodies[i].masses)
else:
velo.append(bodies[i].velocities)
pos.append(bodies[i].trajectory[-1])
mass.append(bodies[i].mass)
n += 1
pos = np.asarray(pos)
velo = np.asarray(velo)
mass = np.asarray(mass)
#return
return pos, velo, mass, n, t, barnes_hut, theta, approximation
def force_calc(mass, positions, D):
......@@ -228,7 +285,7 @@ def force_approx(mass, positions, D, n):
return Force, U
def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave):
def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave, continue_calc, filename = 'none'):
"""
Calculates the dynamics (trajectories) of all bodies
......@@ -252,6 +309,10 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave
every j_th timestep the data is saved
startsave: int
start saving at this timestep
continue_calc: bool
states whether to continue on previous calculation or not
filename:
name of the file of previous calculations
Returns:
--------
......@@ -266,7 +327,11 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave
raise ValueError("Can't combine Barnes-Hut with our approximation. Set at least one to false")
T = round(t_max/h)
positions, velo, mass, n = extract_data(bodies, D)
if continue_calc:
positions, velo, mass, n, t_init, barnes_hut, theta, approximation = extract_file_data(filename)
elif not continue_calc:
positions, velo, mass, n = extract_data(bodies, D)
t_init = 0
if barnes_hut and theta != 0 and not approximation:
force0, U = force_barneshut(positions, mass, theta)
......@@ -283,7 +348,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave
index = -1
if startsave == 0:
index += 1
times.append(0)
times.append(t_init)
E_kin[index] = 0.5*(np.linalg.norm(velo, axis=1)**2)
E_pot[index] = U
all_pos[index] = positions
......@@ -299,7 +364,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave
all_pos[index] = positions
E_kin[index] = 0.5*(np.linalg.norm(velo, axis=1)**2)
E_pot[index] = U
times.append(h*i)
times.append(h*i + t_init)
elif approximation and not barnes_hut:
for i in range(1,T+1):
positions = position(positions, velo, force0, h)
......@@ -311,7 +376,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave
all_pos[index] = positions
E_kin[index] = 0.5*(np.linalg.norm(velo, axis=1)**2)
E_pot[index] = U
times.append(h*i)
times.append(h*i + t_init)
else:
for i in range(1,T+1):
positions = position(positions, velo, force0, h)
......@@ -323,7 +388,7 @@ def dynamics(bodies, D, h, t_max, barnes_hut, theta, approximation, j, startsave
all_pos[index] = positions
E_kin[index] = 0.5*(np.linalg.norm(velo, axis=1)**2)
E_pot[index] = U
times.append(h*i)
times.append(h*i + t_init)
all_pos_corr = np.transpose(np.transpose(all_pos, (1,0,2)) - np.transpose(all_pos, (1,0,2))[0], (1,0,2))
......
......@@ -567,6 +567,17 @@ class Settings:
self.startsave = sets['start_save']
except:
self.startsave = 0
try:
self.continuecalc = sets['continue_calc']
except:
self.continuecalc = False
try:
self.filename = sets['file_name']
except:
self.filename = 'none'
self.continuecalc = 'False'
......
......@@ -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' : 100, # number of asteroids
'N' : 5000, # 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
......@@ -163,7 +163,7 @@ Asteroid_belt = Asteroids(Asteroid_belt_parameters); solar_system.append(Asteroi
# Resonance Astroid belt
Resonance_belt_parameters = {
'N' : 10, # number of asteroids
'N' : 5000, # number of asteroids
'mass' : 0.012, # total mass of all asteroids combined
'R' : 3.1, # orbital radius from the Sun
'theta' : 0,
......
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