Example provided by Deborah Lokhorst. In this example, the four Galilean moons of Jupiter are downloaded from HORIZONS and their orbits are integrated forwards in time. This is a well-known example of a 1:2:4 resonance (also called Laplace resonance) in orbiting bodies. We calculate the resonant arguments see them oscillate with time. We also perform an Fast Fourier Transform (FFT) on the x-position of Io, to look for the period of oscillations caused by the 2:1 resonance between Io and Europa.
Let us first import REBOUND, numpy and matplotlib. We then download the current coordinates for Jupiter and its moons from the NASA HORIZONS database. We work in units of AU, days and solar masses.
In [24]:
import rebound
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
sim = rebound.Simulation()
sim.units = ('AU', 'days', 'Msun')
# We can add Jupiter and four of its moons by name, since REBOUND is linked to the HORIZONS database.
labels = ["Jupiter", "Io", "Europa","Ganymede","Callisto"]
sim.add(labels)
Let us now calculate the mean motions and periods of the inner three moons.
In [2]:
os = sim.calculate_orbits()
print("n_i (in rad/days) = %6.3f, %6.3f, %6.3f" % (os[0].n,os[1].n,os[2].n))
print("P_i (in days) = %6.3f, %6.3f, %6.3f" % (os[0].P,os[1].P,os[2].P))
We can see that the mean motions of each moon are twice that of the moon inner to it and the periods of each moon are half that of the moon inner to it. This means we are close to a 4:2:1 resonance.
Let's move to the center of mass (COM) frame and plot the orbits of the four moons around Jupiter:
In [3]:
sim.move_to_com()
fig = rebound.OrbitPlot(sim, unitlabel="[AU]", color=True, trails=True, periastron=True)
Note that REBOUND automatically plots Jupiter as the central body in this frame, complete with a star symbol (not completely representative of this case, but it'll do).
We can now start integrating the system forward in time. This example uses the symplectic Wisdom-Holman type whfast integrator since no close encounters are expected. The timestep is set to 5% of one of Io's orbits.
In [4]:
sim.integrator = "whfast"
sim.dt = 0.05 * os[0].P # 5% of Io's period
Nout = 100000 # number of points to display
tmax = 80*365.25 # let the simulation run for 80 years
Nmoons = 4
Similar to as was done in the Fourier analysis & resonances example, we set up several arrays to hold values as the simulation runs. This includes the positions of the moons, eccentricities, mean longitudes, and longitude of pericentres.
In [5]:
x = np.zeros((Nmoons,Nout))
ecc = np.zeros((Nmoons,Nout))
longitude = np.zeros((Nmoons,Nout))
varpi = np.zeros((Nmoons,Nout))
times = np.linspace(0.,tmax,Nout)
ps = sim.particles
for i,time in enumerate(times):
sim.integrate(time)
# note we use integrate() with the default exact_finish_time=1, which changes the timestep near
# the outputs to match the output times we want. This is what we want for a Fourier spectrum,
# but technically breaks WHFast's symplectic nature. Not a big deal here.
os = sim.calculate_orbits()
for j in range(Nmoons):
x[j][i] = ps[j+1].x
ecc[j][i] = os[j].e
longitude[j][i] = os[j].l
varpi[j][i] = os[j].Omega + os[j].omega
If we plot the eccentricities as a function of time, one can see that they oscillate significantly for the three inner moons, which are in resonance with each other. Contrasting with these large oscillations, is the smaller oscillation of the outer Galilean moon, Callisto, which is shown for comparison. The three inner moons are in resonance, 1:2:4, but Callisto is not quite in resonance with them, though it is expected to migrate into resonance with them eventually.
Also visible is the gradual change in eccentricity as a function of time: Callisto's mean eccentricity is decreasing and Ganymede's mean eccentricity is increasing. This is a secular change due to the interactions with the inner moons.
In [6]:
fig = plt.figure(figsize=(12,5))
ax = plt.subplot(111)
plt.plot(times,ecc[0],label=labels[1])
plt.plot(times,ecc[1],label=labels[2])
plt.plot(times,ecc[2],label=labels[3])
plt.plot(times,ecc[3],label=labels[4])
ax.set_xlabel("Time (days)")
ax.set_ylabel("Eccentricity")
plt.legend();
We can plot their x-locations as a function of time as well, and observe their relative motions around Jupiter.
In [7]:
fig = plt.figure(figsize=(12,5))
ax = plt.subplot(111)
plt.plot(times,x[0],label=labels[1])
plt.plot(times,x[1],label=labels[2])
plt.plot(times,x[2],label=labels[3])
plt.plot(times,x[3],label=labels[4])
ax.set_xlim(0,0.2*365.25)
ax.set_xlabel("Time (days)")
ax.set_ylabel("x locations (AU)")
ax.tick_params()
plt.legend();
Resonances are identified by looking at the resonant arguments, which are defined as: $$ \theta = (p + q)\lambda_{\rm out} - p \lambda_{\rm in} - q \omega_{\rm out/in}$$ where $\lambda_{\rm out}$ and $\lambda_{\rm in}$ are the mean longitudes of the outer and inner bodies, respectively, and $\omega_{\rm out}$ is the longitude of pericenter of the outer/inner body. The ratio of periods is defined as : $$P_{\rm in}/P_{\rm out} ~= p / (p + q)$$
If the resonant argument, $\theta$, oscillates but is constrained within some range of angles, then there is a resonance between the inner and outer bodies. We call this libration of the angle $\theta$. The trick is to find what the values of q and p are. For our case, we can easily see that there are two 2:1 resonances between the moons, so their resonant arguments would follow the function: $$\theta = 2 \lambda_{\rm out} - \lambda_{\rm in} - \omega_{\rm out}$$ To make the plotting easier, we can borrow this helper function that puts angles into 0 to 360 degrees from another example (Fourier analysis & resonances), and define a new one that puts angles into -180 to 180 degrees.
In [10]:
def zeroTo360(val):
while val < 0:
val += 2*np.pi
while val > 2*np.pi:
val -= 2*np.pi
return (val*180/np.pi)
def min180To180(val):
while val < -np.pi:
val += 2*np.pi
while val > np.pi:
val -= 2*np.pi
return (val*180/np.pi)
# We can calculate theta, the resonant argument of the 1:2 Io-Europa orbital resonance,
# which oscillates about 0 degrees:
theta = [min180To180(2.*longitude[1][i] - longitude[0][i] - varpi[0][i]) for i in range(Nout)]
# There is also a secular resonance argument, corresponding to the difference in the longitude of perihelions:
# This angle oscillates around 180 degs, with a longer period component.
theta_sec = [zeroTo360(-varpi[1][i] + varpi[0][i]) for i in range(Nout)]
fig = plt.figure(figsize=(12,5))
ax = plt.subplot(111)
ax.plot(times,theta)
ax.plot(times,theta_sec) # secular resonance argument
ax.set_xlim([0,20.*365.25])
ax.set_ylim([-180,360.])
ax.set_xlabel("time (days)")
ax.set_ylabel(r"resonant argument $\theta_{2:1}$")
ax.plot([0,100],[180,180],'k--')
ax.plot([0,100],[0,0],'k--')
Out[10]:
Io, Europa and Ganymede are in a Laplace 1:2:4 resonance, which additionally has a longer period libration argument that depends on all three of their mean longitudes, that appears slightly in the other resonant arguments:
In [19]:
thetaL = [zeroTo360(-longitude[0][i] + 3.*longitude[1][i] - 2.*longitude[2][i]) for i in range(Nout)]
fig = plt.figure(figsize=(12,5))
ax = plt.subplot(111)
ax.plot(times,thetaL)
ax.set_ylim([0,360.])
ax.set_xlabel("time (days)")
ax.set_ylabel(r"libration argument $\theta_{2:1}$")
ax.plot([0,200],[180,180],'k--')
Out[19]:
For completeness, let's take a brief look at the Fourier transforms of the x-positions of Io, and see if it has oscillations related to the MMR. We are going to use the scipy Lomb-Scargle periodogram function, which is good for non-uniform time series analysis. Therefore, if we used the IAS15 integrator, which has adaptive timesteps, this function would still work.
In [20]:
from scipy import signal
Npts = 3000
# look for periodicities with periods logarithmically spaced between 0.01 yrs and 100 yrs
logPmin = np.log10(0.001*365.25)
logPmax = np.log10(10.*365.25)
# set up a logspaced array from 0.01 to 100 yrs
Ps = np.logspace(logPmin,logPmax,Npts)
# calculate an array of corresponding angular frequencies
ws = np.asarray([2*np.pi/P for P in Ps])
# calculate the periogram (for Io) (using ws as the values for which to compute it)
periodogram = signal.lombscargle(times,x[0],ws)
In [22]:
fig = plt.figure(figsize=(12,5))
ax = plt.subplot(111)
# Since the computed periodogram is unnormalized, taking the value A**2*N/4,
# we renormalize the results by applying these functions inversely to the output:
ax.set_xscale('log')
ax.set_xlim([10**logPmin,10**logPmax])
ax.set_xlabel("Period (days)")
ax.set_ylabel("Power")
ax.plot(Ps,np.sqrt(4*periodogram/Nout))
Out[22]:
The first spike at about 2 days is caused by the motion of Io around Jupiter in its orbit.
The other spikes, corresponding to oscillations with periods of around 1 year, are caused by the MMR of the moons. The largest spike at ~1.3 years is the from the 1:2 resonance of the two inner moons, Io and Europa.
In [ ]: