Generalized N-Body Solver

This application is designed to solve an n-body gravitational system for two or more bodies. I have designed all of this code around a Simple Example from Towards Data Science. The example I used lacked any sort of class structure and code to animate 3-body solutions. Additionally I've added features such as error inspection and generalized center of mass calculations.

This code has been tested for systems of 2-5 bodies. Expect the exectution time to slow down for more bodies & shorter timesteps.

How to Use

1. Define some galactic body objects

Give each body a name, mass, initial position, initial velocity

alpha_centauri_a = NamedBody('Alpha Centauri A', m=1.1, r0=np.array([-0.5, 0, 0]), v0=np.array([0.01, 0.01, 0]))
alpha_centauri_b = NamedBody('Alpha Centauri B', m=0.907, r0=np.array([0.5, 0, 0]), v0=np.array([-0.05, 0, -0.1]))
body_c = NamedBody('Body C', m=1.0, r0=[0, 1, 0], v0=[0, -0.01, 0])
two_bodies = [alpha_centauri_a, alpha_centauri_b]
three_bodies = two_bodies + [body_c]
input_bodies = two_bodies

2. Initialize number of periods and timesteps
 a. for a single system:

p, step = 50, 0.01  
System = NBodySystem(copy.deepcopy(input_bodies), nd, t_step=step, periods=p)  

 b. For multiple systems:

p, timesteps = 50, [0.01, 0.1]  
Systems = [NBodySystem(copy.deepcopy(input_bodies), nd, t_step=step, periods=p)  for step in timesteps]  

3. Solve each system
The solve loop iterates over each system and applies the following:

• execute() - runs the ODE solve routine and updates the body objects with solutions
• compute_center_of_mass() - computes the system CoM at each timestep from the solutions
• compute_relative_positions() - computes the position of each body relative to the system CoM
• save_solutions( ... ) - save the position and velocity data to a csv file. Provide a directory and unique tag for each solution
  
  

Example (Note this loop functions with 1 or more systems):

save_ = False
for i, system in enumerate(Systems):
	system.execute()
	system.compute_center_of_mass()
	system.compute_relative_positions()
	system.save_solutions(save=save_, dir_='results//two-body', tag='-' + str(len(system.bodies)) + '-' + str(p) + '-' + str(int(1/timesteps[i])))

4. Extract desired data to plot
The solutions are stored within each body within the n-body-system object. Data may be extracted for plotting as so:

position_data = [{body.name: list(body.r_sol.values()) for body in S.bodies} for S in Systems]
velocity_data = [{body.name: list(body.v_sol.values()) for body in S.bodies} for S in Systems]
pos_relative_data = [{body.name: body.r_com for body in S.bodies} for S in Systems]

5. Visualize Solutions - Examples:

Several examples are provided to demonstrate various use-cases

 a. Plot position or velocity of a single system's solution  

The example accepts of dict of {name: data} pairs where the data is the solution for the given body name

EXAMPLE_single_plot(position_data[0])

 b. Animate the position or velocity of a single system

animate_solution(data_pos=position_data[0], data_com=None)

Video of Animated Solution

 c. Visualize position & velocity data for a System on a single figure

EXAMPLE_pos_vel_plot(position_data[0], velocity_data[0], timesteps[0]))

 d. Compare position and velocity visually at two different time incriments

EXAMPLE_compare_timesteps(position_data, velocity_data, timesteps)

Don't forget plt.show(), as it's not included in the examples.