REBOUNDx trouble with satellites

[marshacker]

I'm attempting to get REBOUNDx to apply a direct force to a simple 3-body system. And I'm not understanding something basic.

particle0 is the Sun particle1 is Earth at 1AU particle2 is the Moon at 60-Earth radii from particle1

I'm using the IAS15 integrator with a constant time step dt = 1e-3 years, which is ~1% the lunar orbital period.

---

I'm using the following to tell REBOUNDx to change the Moon's elements computed with respect to the Earth.

mod.params["coordinates"] = reboundx.coordinates["PARTICLE"] ps[2].params["primary"] = 1 (is this right??)

Then I'm setting the tau_a parameter for the Moon and integrate for 100 years.

ps[2].params["tau_a"] = 1000. #some number to get some decay after a century print(ps[2].calculate_orbit(primary=ps[1]).a) #elements before integration sim.integrate(100) print(ps[2].calculate_orbit(primary=ps[1]).a) #elements after integration

But comparing the two prints shows that nothing happens (i.e., the force that gives rise to "tau_a" = 1000 is not applied) unless​ I change the above statement to

ps[1].params["tau_a"] = 1000.

---

It seems that tau_a must be set for Earth to get the Moon's orbit to change.

This doesn't seem right.

Darren Williams

[dtamayo]

Hi Darren, I'm sorry for the long delay. I think this is working correctly, but it is tricky getting the output relative to the right body. Here's how I set things up with roughly a coplanar sun earth moon system:

import rebound
import reboundx
import numpy as np
sim = rebound.Simulation()
sim.G = 4*np.pi**2
sim.add(m=1.)
sim.add(m=3e-6,a=1)
sim.add(m=3e-8,a=60*6.4e3/1.5e8, primary=sim.particles[1])
sim.move_to_com() # Moves to the center of momentum frame
ps = sim.particles

Let's make sure we can read the moon's earth-centric orbit right

o = ps[2].calculate_orbit(primary=ps[1])
a0 = o.a
print(a0, o.e)

This prints 0.00256, 1.8e-15, with the semi major axis of the moon in AU. Now let's set up reboundx:

rebx = reboundx.Extras(sim)
mof = rebx.load_force("modify_orbits_forces")
rebx.add_force(mof)

tmax = 1.e3
ps[2].params["tau_a"] = -tmax*3
mof.params['coordinates'] = reboundx.coordinates["PARTICLE"]
ps[1].params["primary"] = 1

I've made up numbers to integrate for 1000 earth years, with a semi major axis decay timescale 3x longer. Now we run and store the earth-centric semi major axis at every output

Nout = 1000
a = np.zeros(Nout)
times = np.linspace(0.,tmax,Nout)
for i,time in enumerate(times):
    sim.integrate(time)
    orb = ps[2].calculate_orbit(primary=ps[1])
    a[i] = orb.a

and plot

apred = [a0*np.e**(t/ps[2].params["tau_a"]) for t in times]

%matplotlib inline
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(15,5))
ax = plt.subplot(111)
ax.set_yscale('log')
ax.plot(times,a, 'r.',  label='Integration')
ax.plot(times,apred, 'k--',label='Prediction')
ax.set_xlabel("Time", fontsize=24)
ax.set_ylabel("Semimajor axis", fontsize=24)
ax.legend(fontsize=24)

This gives me a semi major axis evolution for the moon that follows the predicted exponential decay. Hope this is helpful (and still relevant!)

[dtamayo]

I'll close this for now for my own checklist, but feel free to reopen it with any follow-up issues/questions. I should be much more responsive after our move this Sunday