###################
# Python (C-like) #
###################
import numpy as np
import sys

# Read in number of timesteps
N_times = int(sys.argv[1])
# Constants
G = 6.67408e-11
m = 42
dt = 1e-3
dt_over_m = dt/m
dt_G_m_m = dt*G*m*m
# Create random initial conditions
N = 100
r = np.random.rand(N, 3)
p = np.random.rand(N, 3)
# Allocations
dr = np.empty(3)

def timestep():
    # Pairwise loop
    for i in range(N):
        ri, pi = r[i], p[i]
        for j in range(i + 1, N):
            rj, pj = r[j], p[j]
            # Gravity
            norm3 = 0
            for dim in range(3):
                dr[dim] = ri[dim] - rj[dim]
                norm3 += dr[dim]**2
            norm3 **= 1.5
            for dim in range(3):
                dp = dt_G_m_m*dr[dim]/norm3
                # Time evolution
                pi[dim] += dp
                pj[dim] -= dp
                ri[dim] += pi[dim]*dt_over_m
                rj[dim] += pj[dim]*dt_over_m

# Time loop
for i in range(N_times):
    timestep()