"""Basic testing to ensure a galaxy evolves unperturbed
in isolation."""

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import savgol_filter

from analysis import utils

f, ax = plt.subplots(1, 1)

time = []
median = []
areas = []

truth = utils.loadData('isolated', 100)
for n in np.linspace(100, 10000, 100):
	
	# Load the data
	data = utils.loadData('isolated', int(n))
	# Find the deviation
	time.append(data['t'])
	deviations = np.abs(np.linalg.norm(data['r_vec'], axis=-1) 
		- np.linalg.norm(truth['r_vec'], axis=-1))
	median.append(np.median(deviations))
	areas.append(np.percentile(deviations, [6, 16, 31, 100-31, 100-16, 100-6]))

# Plot
areas = np.array(areas)
areas = savgol_filter(areas, 11, 3, axis=0) #Smoothing
median = savgol_filter(median, 11, 3)
ax.fill_between(time, areas[:,0], areas[:,-1], alpha=.25, color='darkorange')
ax.fill_between(time, areas[:,1], areas[:,-2], alpha=.25, color='darkorange')
ax.fill_between(time, areas[:,2], areas[:,-3], alpha=.25, color='darkorange')
ax.plot(time, median, color='darkorange')
utils.stylizePlot([ax])
utils.setAxes(ax, x='Time', y='Error in radius', 
	xcoords=(.85,-0.08), ycoords=(-0.1,.75))
ax.set_yscale('log')
ax.set_xlim(0, time[-1])
ax.set_ylim(1E-5, None)

#plt.legend(loc='upper left')
plt.show()