How to use the rebound.get_particles function in rebound

def energy():
        particles = rebound.get_particles()
        com_vx = 0.
        com_vy = 0.
        com_vz = 0.
        if integrator=="wh" or integrator=="mercury":
            mtot = 0.
            for i in xrange(0,N):
                com_vx += particles[i].vx*particles[i].m 
                com_vy += particles[i].vy*particles[i].m 
                com_vz += particles[i].vz*particles[i].m 
                mtot += particles[i].m
            com_vx /= mtot
            com_vy /= mtot
            com_vz /= mtot
        E_kin = 0.
        E_pot = 0.
        for i in xrange(N):

def energy():
        particles = rebound.get_particles()
        com_vx = 0.
        com_vy = 0.
        com_vz = 0.
        if integrator=="wh" or integrator=="mercury":
            mtot = 0.
            for i in xrange(0,N):
                com_vx += particles[i].vx*particles[i].m 
                com_vy += particles[i].vy*particles[i].m 
                com_vz += particles[i].vz*particles[i].m 
                mtot += particles[i].m
            com_vx /= mtot
            com_vy /= mtot
            com_vz /= mtot
        E_kin = 0.
        E_pot = 0.
        for i in xrange(N):

def move_to_heliocentric():
        particles = rebound.get_particles()
        
        particles[0].x  = 0.
        particles[0].y  = 0. 
        particles[0].z  = 0. 
        particles[0].vx = 0. 
        particles[0].vy = 0. 
        particles[0].vz = 0.

def move_to_heliocentric():
        particles = rebound.get_particles()
        
        particles[0].x  = 0.
        particles[0].y  = 0. 
        particles[0].z  = 0. 
        particles[0].vx = 0. 
        particles[0].vy = 0. 
        particles[0].vz = 0.

def move_to_heliocentric():
        particles = rebound.get_particles()
        
        for i in xrange(1,N):
            particles[i].x -= particles[0].x
            particles[i].y -= particles[0].y
            particles[i].z -= particles[0].z
            particles[i].vx -= particles[0].vx
            particles[i].vy -= particles[0].vy
            particles[i].vz -= particles[0].vz
        particles[0].x  = 0.
        particles[0].y  = 0. 
        particles[0].z  = 0. 
        particles[0].vx = 0. 
        particles[0].vy = 0. 
        particles[0].vz = 0.

def energy():
        if integrator=="wh":
            rebound.move_to_center_of_momentum()
        particles = rebound.get_particles()
        E_kin = 0.
        E_pot = 0.
        for i in xrange(N):
            E_kin += 0.5*particles[i].m*(particles[i].vx*particles[i].vx + particles[i].vy*particles[i].vy + particles[i].vz*particles[i].vz)
            for j in xrange(i+1,N):
                dx = particles[i].x-particles[j].x
                dy = particles[i].y-particles[j].y
                dz = particles[i].z-particles[j].z
                r2 = dx*dx + dy*dy + dz*dz
                E_pot -= G*particles[i].m*particles[j].m/np.sqrt(r2)
        if integrator=="wh":
            move_to_heliocentric()
        return E_kin+E_pot