title: N-Body Simulation of Globular Cluster
body: <p>I am trying to code an n-body simulation using the leap frog scheme to simulate a globular cluster however am running into the issue of particles being flung out of the system from (I think) approaching to close to the other particles as well as this leading to a large increase in potential energy. I have tried to write some code to account for collisions, but am unsure what collision radius to use if I am using a softening factor as well.</p>
  <p>The initial positions are randomly generated within a sphere and all have 0 initial velocities.</p>
  <p>If the softening factor is set to 0.1 * Radius of the sphere of generation and the collision radius set to the radius of the sun (leading to no collisions) with N bodies of solar mass I get the following results:</p>
  <p><a href="https://i.stack.imgur.com/QrxPp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/QrxPp.png" alt="paths of bodies over 50 years" /></a>
  <a href="https://i.stack.imgur.com/jECe9.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/jECe9.png" alt="change in energy over 50 years" /></a>
  <a href="https://i.stack.imgur.com/vtRKC.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/vtRKC.png" alt="percentage change of energy over 50 years" /></a></p>
  <p>From the graphs of energy, we can see there a huge spike in potential energy when the bodies originally get close to one another. And over time my energy trails off to 0, it think just due to numerical calculations.</p>
  <p>Is there a way to stop these stars being flung out. I am trying to investigate how long it takes the initial cluster to form an equilibrium state.</p>
  <p>sphere generation:</p>
  <pre><code>sun_mass = 1.989e30
  
  
  N = 100
  mass = sun_mass
  M = np.full([N],mass)
  R = 1e13
  epsilon = 0.1 * R
  collision_radius = 7e8
  V = np.zeros([N,3])
  M = np.full([N],mass)
  P = np.zeros([N,3])
  
  
  t = 50 * 365 * 24 * 60 * 60
  dt = 30 * 24 * 60 * 60
  
  
  
  print(t/dt)
  
  np.random.seed(54321)
  
  for i in range(N):
      
      phi = np.random.uniform(0,2*np.pi)
      costheta = np.random.uniform(-1,1)
      u = np.random.uniform(0,1)
      
      theta = np.arccos( costheta )
      r = R * (u) **(1/3)
      
      
      x = r * np.sin( theta) * np.cos( phi )
      y = r * np.sin( theta) * np.sin( phi )
      z = r * np.cos( theta )
  
      P[i] = (x,y,z)
  </code></pre>
  <p>Programe:</p>
  <pre><code>def programe(position, mass, velocity, softening, time, dt, R, collision_radius):
      
      no_of_time_steps = (round(time/dt))
  
      all_positions = []
      all_velocities = []
      #print(all_positions)
      #print(len(all_positions[0]))
      kinetic_energy = []
      potential_energy = []
      total_energy = []
      
          
      for i in range(no_of_time_steps):
          
          position, mass, velocity = detect_collisions(position, mass, velocity, collision_radius)
          
          all_positions.append(position)
          all_velocities.append(velocity)
      
          'graph'
          plots = np.round(np.linspace(0,no_of_time_steps,num=500))
          for k in range(len(plots)):
              if i == plots[k]:
                  print(&quot;test&quot;)
                  #print(i)
                  graph(position, R, k)
          
          'energies'
          kinetic_energy.append(calc_kinetic_energy(velocity, mass))
          potential_energy.append(calc_potential_energy(position, mass))
          total_energy.append(calc_total_energy(position, velocity, mass))
          
          'leap frog'
          velocity = calc_next_v_half(position, mass, velocity, softening, dt)    
          position = calc_next_position(position, mass, velocity, dt)    
          velocity = calc_next_v_half(position, mass, velocity, softening, dt)
          
  
      
      all_positions = np.array(all_positions)
      graphing(all_positions, time, dt, kinetic_energy, potential_energy, total_energy, no_of_time_steps, R)
      #print(len(mass))
  
      return(all_positions, all_velocities, kinetic_energy, potential_energy, total_energy)
  </code></pre>
  <p>collision detection:</p>
  <pre><code>def indexOf(Array, item):
      for i in range(len(Array)):
          if (Array[i] == item).all():
              return i
          
  def detect_collisions(position, mass, velocity, collision_radius):
      i = 0
      
      newP = position
      newM = mass
      newV = velocity
      #print(len(position), len(newM))
      while i &lt; len(newM):
          if newM[i] == 0:
              i += 1
              continue
          j = i + 1
          while j &lt; len(newP):
             #print(j)
              if newM[j] == 0:
                  j += 1
                  continue
              p1 = position[i]
              p2 = position[j]
              
              if calc_seperation(p1,p2) &lt; collision_radius:
                  index1 = indexOf(position, p1)
                  index2 = indexOf(position, p2)
                  print('collision', index1, index2)
                  newM, newV, newP = handle_collision(newM, newV, newP, [index1, index2])
              j += 1
          i += 1
      
      return(newP, newM, newV)
  
  def handle_collision(M, V, P, indexes):
      if M[indexes[0]] &gt; M[indexes[1]]:
          primary = indexes[0]
          secondary = indexes[1]
      else:
          primary = indexes[1]
          secondary = indexes[0]
  
      primaryMomentum = M[primary] * V[primary]
      secondaryMomentum = M[secondary] * V[secondary]
      newMomentum = primaryMomentum + secondaryMomentum
      
      newMass = M[primary] + M[secondary]
      newVelocity = newMomentum / newMass
      M[primary] = newMass
      V[primary] = newVelocity
      M[secondary] = 0
      V[secondary] = 0
      P[secondary] = 0
  
      newM = np.delete(M, secondary, axis=0)
      newV = np.delete(V, secondary, axis=0)
      newP = np.delete(P, secondary, axis=0)
  
      return (newM, newV, newP)
     
  </code></pre>