<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Simulating Saturn's rings
In this example, we will simulate a small patch of Saturn's rings. The simulation is similar to the C example in examples/shearing_sheet.
We first import REBOUND and numpy, then create an instance of the Simulation class to work with.
Step1:  Next up, setting up several constants. We will be simulating a shearing sheet, a box with shear-periodic boundary conditions. This is a local approximation which makes the approximation that the epicyclic frequency $\Omega$ is the same for all particles. 
We work with a value of $\Omega$ that corresponds to a semi-major axis of $a\sim 130000$ km.
Step2:  Next, we need to let REBOUND know about $\Omega$. Within REBOUND $\Omega$ is used by the integrator SEI, the Symplectic Epicycle Integrator (see Rein and Tremaine 2012).
Step3:  Finally, let us define the surface density of the ring and the particle density.
Step4:  The gravitational constant in SI units is
Step5:  We choose a timestep of 1/1000th of the orbital period.
Step6:  We enable gravitational softening to smear out any potential numerical artefacts at very small scales.
Step7:  Next up, we configure the simulation box. By default REBOUND used no boundary conditions, but here we have shear periodic boundaries and a finite simulation domain, so we need to let REBOUND know about the simulation boxsize (note that it is significantly smaller than $a$, so our local approximation is very good. In this example we'll work in SI units.
Step8:  Because we have shear-periodic boundary conditions, we use ghost boxes to simulate the gravity of neighbouring ring patches. The more ghostboxes we use, the smoother the gravitational force accross the boundary. Here, two layers of ghost boxes in the x and y direction are enough (this is a total of 24 ghost boxes). We don't need ghost boxes in the z direction because a rings is a two dimensional system.
Step9:  We can now setup which REBOUND modules we want to use for our simulation. Besides the SEI integrator and the shear-periodic boundary conditions mentioned above, we select the tree modules for both gravity and collisions. This speeds up the code from $O(N^2)$ to $O(N \log(N))$ for large numbers of particles $N$.
Step10:  When two ring particles collide, they loose energy during their the bounce. We here use a velocity dependent Bridges et. al. coefficient of restitution. It is implemented as a python function (a C implementation would be faster!). We let REBOUND know which function we want to use by setting the coefficient_of_restitution function pointer in the simulation instance.
Step11:  To initialize the particles, we will draw random numbers from a power law distribution.
Step12:  Now we can finally add particles to REBOUND. Note that we initialize particles so that they have initially no velovity relative to the mean shear flow.
Step13:  To see what is going on in our simulation, we create a function to plot the current positions of particles and call it once to visualise the initial conditions.
Step14:  We now integrate for one orbital period $P=2\pi/\Omega$.
Step15:  The integration takes a few seconds, then we can visualise the final particle positions.
<ASSISTANT_TASK:>
Python Code:
import rebound
import numpy as np
sim = rebound.Simulation()
Explanation: Simulating Saturn's rings
In this example, we will simulate a small patch of Saturn's rings. The simulation is similar to the C example in examples/shearing_sheet.
We first import REBOUND and numpy, then create an instance of the Simulation class to work with.
End of explanation
OMEGA = 0.00013143527     # [1/s]
Explanation: Next up, setting up several constants. We will be simulating a shearing sheet, a box with shear-periodic boundary conditions. This is a local approximation which makes the approximation that the epicyclic frequency $\Omega$ is the same for all particles. 
We work with a value of $\Omega$ that corresponds to a semi-major axis of $a\sim 130000$ km.
End of explanation
sim.integrator_sei_OMEGA = OMEGA
Explanation: Next, we need to let REBOUND know about $\Omega$. Within REBOUND $\Omega$ is used by the integrator SEI, the Symplectic Epicycle Integrator (see Rein and Tremaine 2012).
End of explanation
surface_density = 400.    # kg/m^2
particle_density = 400.   # kg/m^3
Explanation: Finally, let us define the surface density of the ring and the particle density.
End of explanation
sim.G = 6.67428e-11       # N m^2 / kg^2
Explanation: The gravitational constant in SI units is
End of explanation
sim.dt = 1e-3*2.*np.pi/OMEGA
Explanation: We choose a timestep of 1/1000th of the orbital period.
End of explanation
sim.softening = 0.2       # [m]
Explanation: We enable gravitational softening to smear out any potential numerical artefacts at very small scales.
End of explanation
boxsize = 200.            # [m]
sim.configure_box(boxsize)
Explanation: Next up, we configure the simulation box. By default REBOUND used no boundary conditions, but here we have shear periodic boundaries and a finite simulation domain, so we need to let REBOUND know about the simulation boxsize (note that it is significantly smaller than $a$, so our local approximation is very good. In this example we'll work in SI units.
End of explanation
sim.configure_ghostboxes(2,2,0)
Explanation: Because we have shear-periodic boundary conditions, we use ghost boxes to simulate the gravity of neighbouring ring patches. The more ghostboxes we use, the smoother the gravitational force accross the boundary. Here, two layers of ghost boxes in the x and y direction are enough (this is a total of 24 ghost boxes). We don't need ghost boxes in the z direction because a rings is a two dimensional system.
End of explanation
sim.integrator = "sei"
sim.boundary   = "shear"
sim.gravity    = "tree"
sim.collision  = "tree"
Explanation: We can now setup which REBOUND modules we want to use for our simulation. Besides the SEI integrator and the shear-periodic boundary conditions mentioned above, we select the tree modules for both gravity and collisions. This speeds up the code from $O(N^2)$ to $O(N \log(N))$ for large numbers of particles $N$.
End of explanation
def cor_bridges(r, v):
        eps = 0.32*pow(abs(v)*100.,-0.234)
        if eps>1.:
            eps=1.
        if eps<0.:
            eps=0.
        return eps
sim.coefficient_of_restitution = cor_bridges
Explanation: When two ring particles collide, they loose energy during their the bounce. We here use a velocity dependent Bridges et. al. coefficient of restitution. It is implemented as a python function (a C implementation would be faster!). We let REBOUND know which function we want to use by setting the coefficient_of_restitution function pointer in the simulation instance.
End of explanation
def powerlaw(slope, min_v, max_v):
    y = np.random.uniform()
    pow_max = pow(max_v, slope+1.)
    pow_min = pow(min_v, slope+1.)
    return pow((pow_max-pow_min)*y + pow_min, 1./(slope+1.))
Explanation: To initialize the particles, we will draw random numbers from a power law distribution.
End of explanation
total_mass = 0.
while total_mass < surface_density*(boxsize**2):
    radius = powerlaw(slope=-3, min_v=1, max_v=4)  # [m]    
    mass = particle_density*4./3.*np.pi*(radius**3)
    x = np.random.uniform(low=-boxsize/2., high=boxsize/2.)
    sim.add(
        m=mass,
        r=radius,
        x=x,
        y=np.random.uniform(low=-boxsize/2., high=boxsize/2.),
        z=np.random.normal(),
        vx = 0.,
        vy = -3./2.*x*OMEGA, 
        vz = 0.)
    total_mass += mass
Explanation: Now we can finally add particles to REBOUND. Note that we initialize particles so that they have initially no velovity relative to the mean shear flow.
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.patches as patches
def plotParticles(sim):
    fig = plt.figure(figsize=(8,8))
    ax = plt.subplot(111,aspect='equal')
    ax.set_ylabel("radial coordinate [m]")
    ax.set_xlabel("azimuthal coordinate [m]")
    ax.set_ylim(-boxsize/2.,boxsize/2.)
    ax.set_xlim(-boxsize/2.,boxsize/2.)
    for i, p in enumerate(sim.particles):
        circ = patches.Circle((p.y, p.x), p.r, facecolor='darkgray', edgecolor='black')
        ax.add_patch(circ)
plotParticles(sim)
Explanation: To see what is going on in our simulation, we create a function to plot the current positions of particles and call it once to visualise the initial conditions.
End of explanation
sim.integrate(2.*np.pi/OMEGA)
Explanation: We now integrate for one orbital period $P=2\pi/\Omega$.
End of explanation
plotParticles(sim)
Explanation: The integration takes a few seconds, then we can visualise the final particle positions.
End of explanation
<END_TASK>