[Python] Does this leapfrog method work for the 3 body problem?

I have been trying to make a leapfrog integration to document the variation of the hamiltonian over time for the 3BP, but I never really grasped how to implement it using the normal half-step method so I tried using a variation but I'm not sure if it's correct.

This is the functions I'm using where the variables p, v, m

p = [array([x,y]), array([x,y]), array([x,y])]
v = [array([x,y]), array([x,y]), array([x,y])]
m = [1, 1, 1]


are numpy arrays:

from numpy import sum
from numpy.linalg import norm
from copy import deepcopy

def H(p, v, m): #hamiltonian function
#sum of kinetic energy for all bodies
T = sum([m*norm(v)**2/2 for i in range(3)])
#sum of potential energy between all bodies
V = -sum([m[0-i]*m[1-i]/norm(p[0-i]-p[1-i]) for i in range(3)])
return T + V

def a(p, n, m): #sum of the acceleration arrays for body n and the two other bodies
return m[n-1]*(p[n-1]-p[n])/(norm(p[n-1]-p[n])**3) + m[n-2]*(p[n-2]-p[n])/(norm(p[n-2]-p[n])**3)

def collision(p): #checks for collisions
for i in range(3):
if norm(p[0-i] - p[1-i]) < 0.1:
return True

#leapfrog
def Leapfrog(P, V, dt, steps, m):
p,v=deepcopy(P),deepcopy(V)
H_L = [H(p, v, m)]

for t in range(steps):
atemp = [a(p, i, m) for i in range(3)] #acceleration at time step i

#calc new values
for i in range(3):
p = p + v*dt + 0.5*atemp*dt**2

for i in range(3):
v = v + 0.5*(atemp + a(p, i, m))*dt #acceleration at timestep i+1


if collision(p):
return H_L
H_L.append(H(p, v, m))
return H_L

Continue reading...