I) $$ force = G \frac{m_1 m_{2} direction}{|direction|^3} $$
II) $$ acceleration = \frac{force}{mass}$$
III) $$ position_{t+1} = position_{t} + \frac{1}{5} acceleration_{t}$$
import numpy as np
import k3d
# Set the value of bodies
bodies_amount = 25
# Create 7-dimensional vector for each body, first 3-elements are positions of body,
# next 3 are accelerations in direction of motion and the last element is body mass
bodies = np.random.random_sample((bodies_amount, 7)).astype(np.float32)
# Set the values of bodies positions and accelerations in a half of unit opened ball
bodies[:, 0:6] -= 0.5
# Decrease the values of acceleration by factor of five hundredth
bodies[:, 3:6] *= 0.05
# Increase the masses by adding zero point five and multiplying the sum by one thousand
bodies[:,6] = (bodies[:,6] + 0.5) * 1000
# Set positions of the first body in the center of coordinate system
# Deprive first body acceleration and raise the mass to million of units value
bodies[0,:] = np.array([0, 0, 0, 0, 0, 0, 1e6])
# Normalize the values of bodies positions to set them on a half of unit closed sphere
bodies[1:, 0:3] = bodies[1:, 0:3] / np.linalg.norm(bodies[1:, 0:3], axis=1).reshape(1, bodies_amount-1).T * 0.5
plot = k3d.plot(antialias=True, camera_auto_fit = False, grid_auto_fit = False)
points = k3d.points(bodies[:, 0:3], color=0, point_size=0.05)
plot += points
plot.display()
Output()
# Set the gravitational constant
G = 6.67E-11
# Create accumulators for trace lines and speeds
lines = []
speeds = []
# of each body
for body_index in range(bodies_amount):
lines.append([])
speeds.append([])
for time_step in range(500):
for first_body_index in range(bodies_amount):
resultant_force = np.zeros(3)
for second_body_index in range(bodies_amount):
# Ignore force calculation for the system of two of the same bodies
if first_body_index == second_body_index:
continue
# Calculate the vector connecting initial point in space of first body with terminal point of second
direction = bodies[second_body_index, 0:3] - bodies[first_body_index, 0:3]
# Calculate force between two bodies from formula I
force = G * bodies[first_body_index, 6] * bodies[second_body_index, 6] * direction
force = force / (np.linalg.norm(direction) ** 3)
# Calculate the sum of forces acting on a first body
resultant_force += force
# Recalculate accelerations of bodies II
bodies[first_body_index, 3:6] = bodies[first_body_index, 3:6] + resultant_force / bodies[first_body_index, 6]
for body_index in range(bodies_amount):
# Recalculate positions from formula III
bodies[body_index, 0:3] = bodies[body_index, 0:3] + bodies[body_index, 3:6] * 0.2
# Append trace lines and the norm of acceleration vector of body in point of space
lines[body_index].append(np.copy(bodies[body_index, 0:3]))
speeds[body_index].append(np.linalg.norm(bodies[body_index, 3:6]))
# Change the positions of point on k3d plot
points.positions = np.copy(bodies[:, 0:3])
# Draw color trace lines of bodies motion
for line, speed in zip(lines, speeds):
plot += k3d.line(np.array(line), attribute=speed, shader="mesh", width=0.0025,
color_range=[0,0.1], color_map=k3d.basic_color_maps.Jet)