How to use the rebound.Particle function in rebound

# This script test the orbit conversion routines of REBOUND.
import sys; sys.path.append('../')
import unittest
import rebound
import random
import math

def almost_equal_wrap_2pi(val1,val2, places):
    diff = val2-val1
    diff2 = diff + 2*math.pi if diff < 0 else diff - 2*math.pi
    return True if min(math.fabs(diff), math.fabs(diff2)) < 10**(-places) else False

class cartesian_to_orbital(unittest.TestCase):
    sun = rebound.Particle( m=1.)
    G = 1.
    N_random_tests = 10 # num of random orbits to test in test_rand_r_to_orb_(f or M)
    
    def test_aew2pi(self):
        places=15
        cases = ((0.,10**(-16), True),
                 (0.,10**(-16), True),
                 (0.,10**(-14), False),
                 (0.,10**(-14), False),
                 (0.1,0.2, False))
        for case in cases:
            self.assertIs(almost_equal_wrap_2pi(case[0], 2*math.pi - case[1], places),case[2], '{}'.format(case))
            self.assertIs(almost_equal_wrap_2pi(2*math.pi - case[1], case[0], places),case[2], '{}'.format(case))
            
    def test_keplers_eq(self):
        '''test Kepler's equation'''

ecc   = np.random.rand()*.1
    inc   = np.deg2rad(np.random.rand()*10.)
    Omega = np.random.rand()*2*np.pi
    omega = np.random.rand()*2*np.pi
    f     = np.random.rand()*2*np.pi
    return a, ecc, inc, Omega, omega, f

# Create a simulation with 'nice' units
sim = rebound.Simulation()
sim.units = ["AU", "MSun", "day"]

# Use IAS15 for the adaptive time-stepping
sim.integrator = 'ias15'

# Add the Sun and Jupiter on a circular orbit
Sun = rebound.Particle(m=1.)
sim.add(Sun)
sim.add(primary=Sun, m=1e-3, a=5., id=1) 
Jup = sim.particles[1]

# Add satellites with random eccentricities, inclinations, and 
# joviancentric semimajor axes the between .1 and 2 Hill radii
rHill = Jup.a * (Jup.m/(3.*Sun.m))**(1./3) 
Ntest = 100
i = 0
while i < Ntest:
    # Draw random orbital elements
    testa,testecc,testinc,testOmega,testomega,testf=get_test_orbit(rHill)
    sim.add(primary=Jup, m=0., a=testa, e=testecc, inc=testinc,
            Omega=testOmega, omega=testomega, f=testf, id=i+2)
    i += 1

def simulation(par):
    saturn_a, saturn_e = par
    sim = rebound.Simulation() 
    sim.integrator = "whfast"
    sim.min_dt = 5.
    sim.dt = 1.
    
    # These parameters are only approximately those of Jupiter and Saturn.
    sun     = rebound.Particle(m=1.)
    sim.add(sun)
    jupiter = sim.add(primary=sun,m=0.000954, a=5.204, M=0.600, omega=0.257, e=0.048)
    saturn  = sim.add(primary=sun,m=0.000285, a=saturn_a, M=0.871, omega=1.616, e=saturn_e)

    sim.move_to_com()
    sim.init_megno()
    # Hide warning messages (WHFast timestep too large)
    with warnings.catch_warnings(record=True) as w: 
        warnings.simplefilter("always")
        sim.integrate(1e3*2.*np.pi)

    return [sim.calculate_megno(),1./(sim.calculate_lyapunov()*2.*np.pi)] # returns MEGNO and Lypunov timescale in years

rebound.set_G(k*k)      # Gravitational constant

# Choose the symplectic WHFast integrator 
rebound.set_integrator("whfast")

# This integrator is not adaptive, so we need to set the timestep
rebound.set_dt(40.)     # 40 days (the time unit depends on the unit of G, see above).

# Setup particles (data taken from NASA Horizons)
# This could also be easily read in from a file.
rebound.particle_add( Particle( m=1.00000597682, x=-4.06428567034226e-3, y=-6.08813756435987e-3, z=-1.66162304225834e-6, vx=+6.69048890636161e-6, vy=-6.33922479583593e-6, vz=-3.13202145590767e-9) )  # Sun
rebound.particle_add( Particle( m=1./1047.355,   x=+3.40546614227466e+0, y=+3.62978190075864e+0, z=+3.42386261766577e-2, vx=-5.59797969310664e-3, vy=+5.51815399480116e-3, vz=-2.66711392865591e-6) )  # Jupiter
rebound.particle_add( Particle( m=1./3501.6,     x=+6.60801554403466e+0, y=+6.38084674585064e+0, z=-1.36145963724542e-1, vx=-4.17354020307064e-3, vy=+3.99723751748116e-3, vz=+1.67206320571441e-5) )  # Saturn
rebound.particle_add( Particle( m=1./22869.,     x=+1.11636331405597e+1, y=+1.60373479057256e+1, z=+3.61783279369958e-1, vx=-3.25884806151064e-3, vy=+2.06438412905916e-3, vz=-2.17699042180559e-5) )  # Uranus
rebound.particle_add( Particle( m=1./19314.,     x=-3.01777243405203e+1, y=+1.91155314998064e+0, z=-1.53887595621042e-1, vx=-2.17471785045538e-4, vy=-3.11361111025884e-3, vz=+3.58344705491441e-5) )  # Neptune
rebound.particle_add( Particle( m=0,             x=-2.13858977531573e+1, y=+3.20719104739886e+1, z=+2.49245689556096e+0, vx=-1.76936577252484e-3, vy=-2.06720938381724e-3, vz=+6.58091931493844e-4) )  # Pluto

# Set the center of momentum to be at the origin
rebound.move_to_center_of_momentum()

# Get the particle data
# Note: this is a pointer and will automatically update as the simulation progresses
particles = rebound.particles_get()
# timestep counter
steps = 0 
# Integrate until t=1e6 (unit of time in this example is days)
while rebound.get_t()<1e6:
    rebound.step()
    steps += 1
    # Print particle positions every 100 timesteps
    if steps%100==0:
        for i in range(rebound.get_N()):