1
0

isolation.py 1.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546
  1. """Basic testing to ensure a galaxy evolves unperturbed
  2. in isolation."""
  3. import matplotlib.pyplot as plt
  4. import numpy as np
  5. from scipy.signal import savgol_filter
  6. from analysis import utils
  7. f, ax = plt.subplots(1, 1)
  8. time = []
  9. median = []
  10. areas = []
  11. truth = utils.loadData('isolated', 100)
  12. for n in np.linspace(100, 10000, 100):
  13. # Load the data
  14. data = utils.loadData('isolated', int(n))
  15. # Find the deviation
  16. time.append(data['t'])
  17. deviations = np.abs(np.linalg.norm(data['r_vec'], axis=-1)
  18. - np.linalg.norm(truth['r_vec'], axis=-1))
  19. median.append(np.median(deviations))
  20. areas.append(np.percentile(deviations, [6, 16, 31, 100-31, 100-16, 100-6]))
  21. # Plot
  22. areas = np.array(areas)
  23. areas = savgol_filter(areas, 11, 3, axis=0) #Smoothing
  24. median = savgol_filter(median, 11, 3)
  25. ax.fill_between(time, areas[:,0], areas[:,-1], alpha=.25, color='darkorange')
  26. ax.fill_between(time, areas[:,1], areas[:,-2], alpha=.25, color='darkorange')
  27. ax.fill_between(time, areas[:,2], areas[:,-3], alpha=.25, color='darkorange')
  28. ax.plot(time, median, color='darkorange')
  29. utils.stylizePlot([ax])
  30. utils.setAxes(ax, x='Time', y='Error in radius',
  31. xcoords=(.85,-0.08), ycoords=(-0.1,.75))
  32. ax.set_yscale('log')
  33. ax.set_xlim(0, time[-1])
  34. ax.set_ylim(1E-5, None)
  35. #plt.legend(loc='upper left')
  36. plt.show()