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:
In [1]:
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
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.
In [2]:
rebx = reboundx.Extras(sim)
rf = rebx.load_force("radiation_forces")
rebx.add_force(rf)
rf.params["c"] = 3.e8
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:
In [3]:
ps["Sun"].params["radiation_source"] = 1
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.
In [4]:
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")
Now we add the grains' physical properties. In order for particles to feel radiation forces, we have to set their beta parameter. $\beta$ is tha ratio of the radiation force to the gravitational force from the star (Burns et al. 1979). One can either set it directly:
In [5]:
ps["p1"].params["beta"] = 0.01
or we can calculate it from more fundamental parameters. REBOUNDx has a convenience function that takes the gravitional 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).
In [6]:
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"]))
Now let's run for 100 years (about 3 Saturn orbits), and look at how the eccentricity varies over a Saturn year:
In [8]:
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)
Out[8]:
In [ ]: