orbit.py 1008 B

123456789101112131415161718192021222324252627282930313233343536373839
  1. """Basic testing to ensure a pair of galaxies follows the correct orbit."""
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. from scipy.signal import savgol_filter
  5. from analysis import utils
  6. f, ax = plt.subplots(1, 1)
  7. time = []
  8. deviations = []
  9. for n in np.linspace(100, 10000, 100):
  10. # Load the data
  11. data = utils.loadData('orbit', int(n))
  12. # Find the deviation
  13. # The analytical orbit is r = 1/(1 + cos(phi))
  14. time.append(data['t'])
  15. phi = np.arctan2(data['r_vec'][0,1], data['r_vec'][0,0])
  16. truthr = 1/(1 + np.cos(phi))
  17. print(truthr, np.linalg.norm(data['r_vec'][0,:2]))
  18. deviation = np.abs(np.linalg.norm(data['r_vec'][0,:2]) - truthr)
  19. deviations.append(deviation)
  20. # Plot
  21. deviations = savgol_filter(deviations, 11, 3)
  22. ax.plot(time, deviations, color='darkorange')
  23. utils.stylizePlot([ax])
  24. utils.setAxes(ax, x='Time', y='Error in radius',
  25. xcoords=(.85,-0.08), ycoords=(-0.1,.75))
  26. ax.set_yscale('log')
  27. ax.set_xlim(0, time[-1])
  28. ax.set_ylim(1E-5, None)
  29. #plt.legend(loc='upper left')
  30. plt.show()