How to use the rebound.Simulation function in rebound

def simulation(par):
    integrator, run, trial = par
    sim = rebound.Simulation()
    k = 0.01720209895    
    Gfac = 1./k
    sim.dt = dt
    if integrator == "whfast-nocor":
        integrator = "whfast"
    else:
        sim.ri_whfast.corrector = 11
    sim.integrator = integrator 
    sim.ri_whfast.safe_mode = 0

    massfac = 1.
    sim.add(m=1.00000597682, x=-4.06428567034226e-3, y=-6.08813756435987e-3, z=-1.66162304225834e-6,      vx=+6.69048890636161e-6*Gfac, vy=-6.33922479583593e-6*Gfac, vz=-3.13202145590767e-9*Gfac)   # Sun
    sim.add(m=massfac/1407.355,   x=+3.40546614227466e+0, y=+3.62978190075864e+0, z=+3.42386261766577e-2, vx=-5.59797969310664e-3*Gfac, vy=+5.51815399480116e-3*Gfac, vz=-2.66711392865591e-6*Gfac)   # Jupiter
    sim.add(m=massfac/3501.6,     x=+6.60801554403466e+0, y=+6.38084674585064e+0, z=-1.36145963724542e-1, vx=-4.17354020307064e-3*Gfac, vy=+3.99723751748116e-3*Gfac, vz=+1.67206320571441e-5*Gfac)   # Saturn
    sim.add(m=massfac/22869.,     x=+1.11636331405597e+1, y=+1.60373479057256e+1, z=+3.61783279369958e-1, vx=-3.25884806151064e-3*Gfac, vy=+2.06438412905916e-3*Gfac, vz=-2.17699042180559e-5*Gfac)   # Uranus
    sim.add(m=massfac/19314.,     x=-3.01777243405203e+1, y=+1.91155314998064e+0, z=-1.53887595621042e-1, vx=-2.17471785045538e-4*Gfac, vy=-3.11361111025884e-3*Gfac, vz=+3.58344705491441e-5*Gfac)   # Neptune

import os.path
import os
filename = "cache.bin"

if 'TRAVIS' in os.environ:
    # Shorter built time
    solar_system_objects = ["Sun", "Mercury"]
else:
    # More planets
    solar_system_objects = ["Sun", "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "C/2014 Q2"]

if os.path.isfile(filename):
    # Try to load simulation from file
    sim = rebound.Simulation(filename)
else: 
    sim = rebound.Simulation()
    # Get data from NASA Horizons
    try:
        sim.add(solar_system_objects)
    except socket.error:
        print("A socket error occured. Maybe Horizons is down?")
        sys.exit(0) # we ignore the error and exit

    sim.move_to_com()
    # Configure simulation
    sim.integrator = "whfast"
    sim.set_dt = 0.01
    # Let's save it for next time
    # Note: sim.save() only saves the particle data, not the integrator settings, etc.
    sim.save(filename)  

sim.status()

from compute_satellite_orbits import compute_satellite_orbits as cso
import os

# Function to draw random orbital elements within fixed intervals
def get_test_orbit(rHill):
    # semimajor axis is between .1 and 2 Hill radii
    a     = rHill * (np.random.rand()*(1-.1) + .1)
    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

# Import the rebound module
import rebound
import os

# Create a REBOUND simulation
sim = rebound.Simulation()

# Set variables (defaults are G=1, t=0, dt=0.01)
k = 0.01720209895       # Gaussian constant 
sim.G = k*k         # Gravitational constant

# Setup particles (data taken from NASA Horizons)
# This could also be easily read in from a file.
sim.add( 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
sim.add( 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
sim.add( 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
sim.add( 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
sim.add( 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
sim.add( 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
sim.move_to_com()

# Import the rebound module 
import rebound

# Add particles
# We work in units where G=1.  
sim = rebound.Simulation()
sim.add(m=1. )                  # Test particle
sim.add(m=1e-3,x=1.,vy=1. )     # Planet

# Move particles so that the center of mass is (and stays) at the origin  
sim.move_to_com()

# You can provide a function, written in python to REBOUND.
# This function gets called every time the forces are evaluated.
# Simple add any any additional (non-gravitational) forces to the 
# particle accelerations. Here, we add a simple drag force. This 
# will make the planet spiral into the star.
ps = sim.particles
def dragforce(reb_sim):
    dragcoefficient = 1e-2
    for p in ps:
        p.ax += -dragcoefficient * p.vx

ms = 5, markeredgewidth = 2)
            pt_zy[bi], = axzy.plot(b.z_hr[ti], b.y_hr[ti], 'o', color = 'lightgrey',
                                   alpha = 1, markeredgecolor = b.color,
                                   zorder = 99, clip_on = False,
                                   ms = 5, markeredgewidth = 2)

        # Plot the orbits. We will run rebound
        # backwards in time to get the orbital
        # trails. This avoids having to deal with
        # the case when the occultation is at the
        # beginning of the timeseries!
        if not self.settings.quiet:
            print("Computing orbits for plotting...")
        per = np.max([b.per for b in self.bodies])
        tlong = np.arange(0, trail_dur, self.settings.timestep)
        sim = rebound.Simulation()
        sim.G = GEARTH
        for b in self.bodies:
            sim.add(m = b._m, x = b.x_hr[tf], y = b.y_hr[tf],
                    z = b.z_hr[tf], vx = -b.vx_hr[tf],
                    vy = -b.vy_hr[tf], vz = -b.vz_hr[tf])
        x = np.zeros((len(self.bodies), len(tlong)))
        y = np.zeros((len(self.bodies), len(tlong)))
        z = np.zeros((len(self.bodies), len(tlong)))
        for i, time in enumerate(tlong):
            sim.integrate(time, exact_finish_time = 1)
            for q in range(len(self.bodies)):
                x[q,i] = sim.particles[q].x
                y[q,i] = sim.particles[q].y
                z[q,i] = sim.particles[q].z
        maxx = maxy = maxz = 0.
        for j, b in enumerate(self.bodies):

def generate(objects, Ndays=None, Noutputs=None):
    # create rebound simulation 
    # for object parameters see: 
    # https://rebound.readthedocs.io/en/latest/_modules/rebound/particle.html
    sim = rebound.Simulation()
    sim.units = ('day', 'AU', 'Msun')
    for i in range(len(objects)):
        sim.add( **objects[i] ) 
    sim.move_to_com() 

    if Ndays and Noutputs:
        return integrate(sim, objects, Ndays, Noutputs)
    else:
        return sim

def plot_top_down_view(params_median, params_star, a=None, timestep=None, scaling=30., colors=sns.color_palette('deep'), linewidth=2, plot_arrow=False):
    
    sim = rebound.Simulation()
    sim.add(m=1)
    
    for i, companion in enumerate(config.BASEMENT.settings['companions_all']):
        if (i==0) and (timestep is None): 
            timestep = params_median[companion+'_epoch'] #calculate it for the timestep where the first companion is in transit
        first_epoch = get_first_epoch(timestep, params_median[companion+'_epoch'], params_median[companion+'_period'])
        phase = calc_phase(timestep, params_median[companion+'_period'], first_epoch)
        ecc = params_median[companion+'_f_s']**2 + params_median[companion+'_f_c']**2
        w = np.arccos( params_median[companion+'_f_c'] / np.sqrt(ecc) ) #in rad
        inc = params_median[companion+'_incl']/180.*np.pi
        if a is None:
            a1 = params_star['R_star'] / params_median[companion+'_radius_1'] #in Rsun 
            a1 *= 0.004650467260962157 #in AU
        else:
            a1 = a[i]
#        print(a, inc, ecc, w, phase*2*np.pi)

# Import the rebound module
import rebound
import os

filename = "simulationarchive.bin"
try:
    sim = rebound.Simulation(filename)
    print("Restarting from simulation archive. Last snapshot found at t=%.1f"%sim.t)
except:
    print("Cannot load SimulationArchive. Creating new simulation.")
    sim = rebound.Simulation()
    sim.add(m=1) # star
    sim.add(m=1e-3, a=1, e=0.01) # planet 1
    sim.add(m=1e-3, a=2.5, e=0.01) # planet 2
    sim.integrator = "whfast"
    sim.dt = 3.1415*2.*6./365.25 # 6 days in units where G=1
    sim.move_to_com()

sim.automateSimulationArchive(filename, interval=2.*3.1415*1e5)

# Run a very long simulation.
# Interrupted at any time and then run script again to restart.
sim.integrate(2.*3.1415*1e10) # 10 Gyr