eccentricity.py 1.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445
  1. """Plot the eccentricity of the orbits of the particles in an encounter"""
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. from analysis import utils
  5. from analysis.segmentation import segmentEncounter
  6. # Segment the encounter
  7. data = utils.loadData('200mass', 15000)
  8. (bridgeMask, stolenMask, orbittingMask, tailMask), _, _ = segmentEncounter(data)
  9. # Avoid central masses
  10. orbittingMask, stolenMask = orbittingMask[1:], stolenMask[1:]
  11. # Calculate the eccentricities
  12. tailE = utils.calculateEccentricity(1.0, data['r_vec'][0], data['v_vec'][0],
  13. data['r_vec'][tailMask], data['v_vec'][tailMask])
  14. stolenE = utils.calculateEccentricity(2.0, data['r_vec'][1], data['v_vec'][1],
  15. data['r_vec'][stolenMask], data['v_vec'][stolenMask])
  16. orbittingE = utils.calculateEccentricity(1.0, data['r_vec'][0], data['v_vec'][0],
  17. data['r_vec'][orbittingMask], data['v_vec'][orbittingMask])
  18. # Plotting
  19. f, ax = plt.subplots(1, 1)
  20. bins = np.linspace(0, 10, 301)
  21. style = {'linestyle':'solid', 'normed':True, 'color':'black'}
  22. ax.hist(tailE, bins=bins, histtype='step', label='tail', **style)
  23. ax.hist(stolenE, bins=bins, histtype='step', label='stolen', **style)
  24. ax.hist(orbittingE, bins=bins, histtype='step', label='orbitting', **style)
  25. # Plotting styling
  26. utils.setSize(ax, x=(0,2))
  27. utils.setAxes(ax, x='Eccentricity', y='Density of particles',
  28. xcoords=(.85,-0.08), ycoords=(-0.1,.65))
  29. ax.legend(loc='upper right')
  30. utils.stylizePlot([ax])
  31. ax.axvline(x=1, c='black')
  32. plt.show()
  33. # Summarize the results
  34. print('Not shown, fraction of particles in tail with e>2.0:',
  35. np.sum(tailE > 2)/len(tailE))
  36. print('Tail, higher than e>1:', np.sum(tailE > 1)/len(tailE))
  37. print('Stolen, higher than e>1:', np.sum(stolenE > 1)/len(stolenE))
  38. print('Orbitting, higher than e>1:', np.sum(orbittingE > 1)/len(orbittingE))