12345678910111213141516171819202122232425262728293031323334353637383940414243444546 |
- """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()
-
|