<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Radiation Forces on Circumplanetary Dust
This example shows how to integrate circumplanetary dust particles under the action of radiation forces.  We use Saturn's Phoebe ring as an example, a distant ring of debris. 
We have to make sure we add all quantities in the same units.  Here we choose to use SI units.  We begin by adding the Sun and Saturn, and use Saturn's orbital plane as the reference plane
Step1:  Now let's set up REBOUNDx and add radiation_forces.  We also have to set the speed of light in the units we want to use.
Step2:  By default, the radiation_forces effect assumes the particle at index 0 is the source of the radiation.  If you'd like to use a different one, or it's possible that the radiation source might move to a different index (e.g. with a custom merger routine), you can add a radiation_source flag to the appropriate particle like this
Step3:  Here we show how to add two dust grains to the simulation in different ways.  Let's first initialize their orbits.  In both cases we use the orbital elements of Saturn's irregular satellite Phoebe, which the dust grains will inherit upon release (Tamayo et al. 2011).  Since the dust grains don't interact with one another, putting them on top of each other is OK.
Step4:  Now we add the grains' physical properties.  In order for particles to feel radiation forces, we have to set their beta parameter.  $\beta$ is the ratio of the radiation force to the gravitational force from the star (Burns et al. 1979).  One can either set it directly
Step5:  or we can calculate it from more fundamental parameters.  REBOUNDx has a convenience function that takes the gravitational constant, speed of light, radiation source's mass and luminosity, and then the grain's physical radius, bulk density, and radiation pressure coefficient Q_pr (Burns et al. 1979, equals 1 in the limit that the grain size is >> the radiation's wavelength).
Step6:  Now let's run for 100 years (about 3 Saturn orbits), and look at how the eccentricity varies over a Saturn year
<ASSISTANT_TASK:>
Python Code:
import rebound
import reboundx
import numpy as np
sim = rebound.Simulation()
sim.G = 6.674e-11 # SI units
sim.dt = 1.e4 # Initial timestep in sec.
sim.N_active = 2 # Make it so dust particles don't interact with one another gravitationally
sim.add(m=1.99e30, hash="Sun") # add Sun with mass in kg
sim.add(m=5.68e26, a=1.43e12, e=0.056, pomega = 0., f=0., hash="Saturn") # Add Saturn at pericenter
ps = sim.particles
Explanation: Radiation Forces on Circumplanetary Dust
This example shows how to integrate circumplanetary dust particles under the action of radiation forces.  We use Saturn's Phoebe ring as an example, a distant ring of debris. 
We have to make sure we add all quantities in the same units.  Here we choose to use SI units.  We begin by adding the Sun and Saturn, and use Saturn's orbital plane as the reference plane:
End of explanation
rebx = reboundx.Extras(sim)
rf = rebx.load_force("radiation_forces")
rebx.add_force(rf)
rf.params["c"] = 3.e8
Explanation: Now let's set up REBOUNDx and add radiation_forces.  We also have to set the speed of light in the units we want to use.
End of explanation
ps["Sun"].params["radiation_source"] = 1
Explanation: By default, the radiation_forces effect assumes the particle at index 0 is the source of the radiation.  If you'd like to use a different one, or it's possible that the radiation source might move to a different index (e.g. with a custom merger routine), you can add a radiation_source flag to the appropriate particle like this:
End of explanation
a = 1.3e10 # in meters
e = 0.16
inc = 175*np.pi/180.
Omega = 0. # longitude of node
omega = 0. # argument of pericenter
f = 0. # true anomaly
# Add two dust grains with the same orbit
sim.add(primary=ps["Saturn"], a=a, e=e, inc=inc, Omega=Omega, omega=omega, f=f, hash="p1")
sim.add(primary=ps["Saturn"], a=a, e=e, inc=inc, Omega=Omega, omega=omega, f=f, hash="p2")
Explanation: Here we show how to add two dust grains to the simulation in different ways.  Let's first initialize their orbits.  In both cases we use the orbital elements of Saturn's irregular satellite Phoebe, which the dust grains will inherit upon release (Tamayo et al. 2011).  Since the dust grains don't interact with one another, putting them on top of each other is OK.
End of explanation
ps["p1"].params["beta"] = 0.01
Explanation: Now we add the grains' physical properties.  In order for particles to feel radiation forces, we have to set their beta parameter.  $\beta$ is the ratio of the radiation force to the gravitational force from the star (Burns et al. 1979).  One can either set it directly:
End of explanation
grain_radius = 1.e-5 # grain radius in m
density = 1000. # kg/m^3 = 1g/cc
Q_pr = 1.
luminosity = 3.85e26 # Watts
ps["p2"].params["beta"] = rebx.rad_calc_beta(sim.G, rf.params["c"], ps[0].m, luminosity, grain_radius, density, Q_pr)
print("Particle 2's beta parameter = {0}".format(ps["p2"].params["beta"]))
Explanation: or we can calculate it from more fundamental parameters.  REBOUNDx has a convenience function that takes the gravitational constant, speed of light, radiation source's mass and luminosity, and then the grain's physical radius, bulk density, and radiation pressure coefficient Q_pr (Burns et al. 1979, equals 1 in the limit that the grain size is >> the radiation's wavelength).
End of explanation
yr = 365*24*3600 # s
Noutput = 1000
times = np.linspace(0,100.*yr, Noutput)
e1, e2 = np.zeros(Noutput), np.zeros(Noutput)
sim.move_to_com() # move to center of mass frame first
for i, time in enumerate(times):
    sim.integrate(time)
    e1[i] = ps["p1"].calculate_orbit(primary=ps["Saturn"]).e
    e2[i] = ps["p2"].calculate_orbit(primary=ps["Saturn"]).e
   
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(15,5))
ax.plot(times/yr, e1, label=r"$\beta$={0:.1e}".format(ps["p1"].params["beta"]))
ax.plot(times/yr, e2, label=r"$\beta$={0:.1e}".format(ps["p2"].params["beta"]))
ax.set_xlabel('Time (yrs)', fontsize=24)
ax.set_ylabel('Eccentricity', fontsize=24)
plt.legend(fontsize=24)
Explanation: Now let's run for 100 years (about 3 Saturn orbits), and look at how the eccentricity varies over a Saturn year:
End of explanation
<END_TASK>