Commit 4e4d881b authored by Olaf's avatar Olaf

Add error to energy plot

parent 946518ae
......@@ -2,6 +2,7 @@ import pickle
import numpy as np
import matplotlib.pyplot as plt
from objects import *
from errorest import *
......@@ -121,6 +122,7 @@ def plot_energy(bodies, settings):
settings: objects.Settings class object
Settings for the interacting system
"""
error_alpha = 0.2
t = bodies[0].times
E_kin = np.zeros((len(t)))
E_pot = np.zeros((len(t)))
......@@ -135,11 +137,22 @@ def plot_energy(bodies, settings):
else:
E_kin += np.asarray(bodies[i].E_kin)
E_pot += np.asarray(bodies[i].E_pot)
E_kin_std, E_kin_mean = error_combine(E_kin)
E_pot_std, E_pot_mean = error_combine(E_pot)
E_std,E_mean = error_combine(E_kin + E_pot)
fig = plt.figure()
plt.plot(t, E_kin, color = 'yellow', label = 'Kinetic energy')
plt.fill_between(t, E_kin - E_kin_std, E_kin + E_kin_std, alpha = error_alpha, color = 'yellow')
plt.plot(t, E_pot, color = 'blue', label = 'Potential energy')
plt.fill_between(t, E_pot - E_pot_std, E_pot + E_pot_std, alpha = error_alpha, color = 'blue')
plt.plot(t, E_kin + E_pot, color = 'green', label = 'Total energy')
plt.fill_between(t, E_kin + E_pot - E_std, E_kin + E_pot + E_std, alpha = error_alpha, color = 'green')
plt.xlabel('time $[y_J]$', fontsize = settings.plot.fontsize_title)
plt.ylabel('Energy $[\\frac{Au^2}{Y^2}]$', fontsize = settings.plot.fontsize_title)
plt.title('Energy')
......
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